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Energy  storage  components  improve  the  energy  efficiency  of  systems  by  reducing  the  mismatch 
between  supply  and  demand.  For  this  purpose,  phase-change  materials  are  particularly  attractive  since 
they  provide  a  high-energy  storage  density  at  a  constant  temperature  which  corresponds  to  the  phase 
transition  temperature  of  the  material.  Nevertheless,  the  incorporation  of  phase-change  materials 
(PCMs)  in  a  particular  application  calls  for  an  analysis  that  will  enable  the  researcher  to  optimize 
performances  of  systems.  Due  to  the  non-linear  nature  of  the  problem,  numerical  analysis  is  generally 
required  to  obtain  appropriate  solutions  for  the  thermal  behavior  of  systems.  Therefore,  a  large  amount 
of  research  has  been  carried  out  on  PCMs  behavior  predictions.  The  review  will  present  models  based  on 
the  first  law  and  on  the  second  law  of  thermodynamics.  It  shows  selected  results  for  several 
configurations,  from  numerous  authors  so  as  to  enable  one  to  start  his/her  research  with  an  exhaustive 
overview  of  the  subject.  This  overview  stresses  the  need  to  match  experimental  investigations  with 
recent  numerical  analyses  since  in  recent  years,  models  mostly  rely  on  other  models  in  their  validation 
stages. 
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1.  Introduction 

The  ever  increasing  level  of  greenhouse  gas  emissions 
combined  with  the  overall  rise  in  fuel  prices  (although  fluctuations 
occur)  are  the  main  reasons  behind  efforts  devoted  to  improve  the 
use  of  various  sources  of  energy.  Economists,  scientists,  and 
engineers  throughout  the  world  are  in  search  of  (1)  strategies  to 
reduce  the  demand;  (2)  methods  to  ensure  the  security  of  the 
supplies;  (3)  technologies  to  increase  the  energy  efficiency  of 
power  systems;  (4)  new  and  renewable  sources  of  energy  to 
replace  the  limited  and  harmful  fossil  fuels. 

One  of  the  options  to  improve  energy  efficiency  is  to  develop 
energy  storage  devices  and  systems  in  order  to  reduce  the 
mismatch  between  supply  and  demand.  Such  devices  and  systems 
also  improve  the  performance  and  reliability  by  reducing  peak 
loads  and  allowing  systems  to  work  within  an  optimal  range. 
Thus,  they  play  a  preponderant  role  in  conserving  energy.  The 
different  forms  of  energy  that  can  be  stored  are  mechanical, 
electrical,  and  thermal.  Here,  mechanical  (gravitational,  com¬ 
pressed  air,  flywheels)  and  electrical  (batteries)  storages  are  not 
considered  while  thermal  energy  storage  is  discussed  in  the 
context  of  latent  heat  (sensible  heat  and  thermochemical  heat  are 
not  considered). 

Latent  heat  storage  is  based  on  the  capture/release  of  energy 
when  a  material  undergoes  a  phase  change  from  solid  to  liquid, 
liquid  to  gas,  or  vice  versa.  Latent  heat  storage  is  particularly 
attractive  since  it  provides  a  high-energy  storage  density  and  has 
the  capacity  to  store  energy  at  a  constant  temperature  -  or  over  a 
limited  range  of  temperature  variation  -  which  is  the  temperature 
that  corresponds  to  the  phase  transition  temperature  of  the 
material.  For  instance,  it  takes  80  times  as  much  energy  to  melt  a 
given  mass  of  water  ( ice)  than  to  raise  the  same  amount  of  water  by 
1  °C.  Table  1  provides  a  typical  comparison  between  properties  of 
different  thermal  storage  materials  used  at  room  temperature.  For 
the  interested  reader,  excellent  global  reviews  that  pertain  to 
phase-change  materials  and  their  various  applications  were 
proposed  by  Farid  et  al.  [1],  Sharma  et  al.  [2],  Zhang  et  al.  [3], 
Regin  et  al.  [4],  Tyagi  and  Buddhi  [5],  Mondal  [6],  Sethi  and  Sharma 
[7],  and  especially  the  recent  one  by  Verma  et  al.  [8]. 

Nevertheless,  the  incorporation  of  phase-change  materials 
(PCMs)  in  a  particular  application  calls  for  an  analysis  that  will 
enable  the  researcher  to  determine  whether  or  not  PCMs  will 
improve  performances  sufficiently  to  justify  extra  costs  for 
additional  systems  and/or  controls  needed.  Mathematical  model¬ 
ing  of  latent  heat  energy  storage  materials  and/or  systems  is 
needed  for  optimal  design  and  material  selection.  Therefore,  a  large 
amount  of  research  has  been  carried  out  on  PCMs  behavior 
predictions  whether  they  are  considered  separately  or  within 
specific  systems. 


Table  1  shows  the  main  differences,  on  average  basis,  between 
sensible  heat  storage  materials  and  latent  heat  storage  classes  of 
materials.  The  table  indicates  that  generally,  PCMs  require  much 
less  mass  or  volume  to  store  the  same  amount  of  energy  at  a  more 
or  less  constant  temperature. 

This  paper  is  building  on  previous  reviews  [1-8]  to  update  the 
available  references  that  pertain  to  mathematical  modeling  and 
simulation  of  thermal  energy  storage  with  phase-change  materi¬ 
als.  First,  it  presents  the  fundamental  mathematical  description  of 
the  phenomenon,  the  Stephan  problem.  Then,  it  provides  basic 
mathematical  descriptions  used  as  basis  for  numerical  modeling 
using  either  first  or  second  law  approaches  and  fixed  or  adaptative 
meshes.  The  next  section,  considered  by  the  authors  as  the  major 
contribution,  is  a  model  collection  of  most  recent  works  published 
on  the  subject.  This  survey  is  organized  according  to  the  problem 
geometry  (Cartesian,  spherical,  and  cylindrical)  and  specific 
configurations  or  applications  (packed  beds,  finned  surfaces, 
porous  and  fibrous  materials,  slurries).  A  synthesis  is  provided 
at  the  end  and  several  recommendations  are  formulated. 

2.  The  Stephan  problem 

Phase  transition  of  a  material  is  described  by  a  particular  kind  of 
boundary  value  problems  for  partial  differential  equations,  where 
phase  boundary  can  move  with  time.  This  question  has  been  first 
studied  by  Clapeyron  and  Lame  in  1831  when  analyzing  the 
formation  of  the  Earth's  crust  when  it  cooled.  In  that  case,  the 
problem  was  simplified  from  a  spherical  geometry  to  a  one¬ 
dimensional  semi-infinite  slab  [9].  This  solution  was  found 
independently  by  Franz  Neumann,  who  introduced  it  in  his 
lectures  notes  of  1835-1840  [10].  Nevertheless,  this  type  of 
problems  is  named  after  Jozef  Stefan,  the  Slovene  physicist  who 
introduced  the  general  class  of  such  problems  in  1889  [11]  in 
relation  to  problems  of  ice  formation.  Existence  of  a  solution  was 
proved  by  Evans  in  1951  [12],  while  the  uniqueness  was  proved  by 
Douglas  in  1957  [13]. 

Very  few  analytical  solutions  are  available  in  closed  form.  They 
are  mainly  for  the  one-dimensional  cases  of  an  infinite  or  semi¬ 
infinite  region  with  simple  initial  and  boundary  conditions  and 
constant  thermal  properties.  Under  these  conditions,  these  exact 
solutions  usually  take  the  form  of  functions  of  the  single  variable  x/ 
t1/2  and  are  known  as  similarity  solutions  [14,15].  A  collection  of 
similarity  solutions  and  references  is  to  be  found  in  [16,17]. 

3.  Numerical  solution 

The  problem  of  predicting  the  behavior  of  phase-change 
systems  is  difficult  due  to  its  inherent  non-linear  nature  at  moving 
interfaces,  for  which  displacement  rate  is  controlled  by  the  latent 


Table  1 

Common  heat  storage  materials. 


Property 

Materials 

Rock 

Water 

Organic  PCM 

Inorganic  PCM 

Density  [kg/m3] 

2240 

1000 

800 

1600 

Specific  heat  [kj/kg  K] 

1.0 

4.2 

2.0 

2.0 

Latent  heat  [kj/kg] 

- 

190 

230 

Latent  heat  [MJ/m3] 

- 

- 

152 

368 

Storage  mass  for  109J,  avg  [kg] 

67,000 

16,000 

5300 

4350 

(At=15K) 

(At= 15  K) 

Storage  volume  for  109J,  avg  [m3] 

30 

16 

6.6 

2.7 

(At=15K) 

(At= 15  K) 

Relative  storage  mass 

15 

4 

1.25 

1.0 

(At=15I<) 

(At= 15  K) 

Relative  storage  volume 

11 

6 

2.5 

1.0 

(At=15K) 

(At= 15  K) 
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heat  lost  or  absorbed  at  the  boundary.  The  following  equation, 
known  as  the  Stephan  condition,  describes  this  process: 


where  X  is  the  latent  heat  of  fusion,  p  is  the  density  (it  is  not 
specified  if  it  is  solid  or  liquid),  s(t)  is  the  surface  position,  k  is  the 
thermal  conductivity,  t  is  time,  and  T  is  the  temperature.  Indexes  s 
and  /  refers  to  solid  and  liquid  phases. 

In  this  situation,  the  position  and  the  velocity  of  boundaries  are 
not  known  a  priori.  In  addition,  since  two  phases  possess  different 
physical  properties,  this  could  create  the  in  numerical  model  non¬ 
physical  discontinuities  which  need  to  be  addressed. 

3.1.  Fixed  grid 

By  introducing  an  enthalpy  method,  the  phase-change  problem 
becomes  much  simpler  since  the  governing  equation  is  the  same 
for  the  two  phases;  interface  conditions  are  automatically  achieved 
and  create  a  mushy  zone  between  the  two  phases.  This  zone  avoids 
sharp  discontinuities  that  may  create  some  numerical  instabilities. 
In  consequence,  the  thickness  and  the  quality  of  the  discretization 
of  this  mushy  zone  are  critical  to  the  model  performance.  The 
enthalpy  method  can  deal  with  both  mushy  and  isothermal  phase- 
change  problems  but  the  temperature  at  a  typical  grid  point  may 
oscillate  with  time  [18].  This  method  has  been  successfully  applied 
to  various  phase-change  problems  [19-22].  Hunter  in  1989  [23] 
and  Amdjadi  in  1990  [24]  confirmed  that  the  enthalpy  method  is 
the  most  suitable  for  typical  applications  under  the  restriction  that 
there  is  no  alteration  to  the  numerical  scheme  at  the  interface. 

Enthalpy  function  h,  defined  as  a  function  of  temperature,  is  given 
by  Voller  [25].  For  a  phase-change  process,  energy  conservation  can 
be  expressed  in  terms  of  total  volumetric  enthalpy  and  temperature 
for  constant  thermophysical  properties,  as  follows  (from  [2]): 

^  =  V[MVT)],  (2) 

where  H  is  the  total  volumetric  enthalpy,  which  is  the  sum  of 
sensible  and  latent  heats: 

H(T)  =  h(T)  +  p,f(T)X  (3) 

and  where 


capacity  method  [27].  In  these  models,  the  enthalpy-based  energy 
equation  is  converted  into  a  non-linear  equation  with  a  single 
variable.  This  approach  is  called  the  temperature  transforming 
model  (TTM).  This  type  of  simulation  has  proved  accurate,  simple 
and  efficient  [27]. 

In  TTM,  general  continuity  and  momentum  equations  for  fluid 
problems  are  used.  But  unlike  in  the  enthalpy  formalism,  the 
energy  equation  uses  a  temperature  dependant  source  term.  The 
main  difficulty  of  this  technique  is  to  develop  a  method  to  keep 
velocity  at  zero  in  the  solid  phase.  The  simplest  method  is  the 
switch-off  method  (SOM)  or  its  smoothed  version  the  ramped 
switch-off  method  (RSOM)  [20,28,29].  These  techniques  directly  set 
velocities  in  the  solid  phase  to  zero  by  setting  coefficients  of 
momentum  and  velocity-correction  equations  to  values  that 
effectively  suppress  any  motion.  Although  these  methods  are 
commonly  used  in  phase-change  problems,  Ma  and  Zhang  [18] 
have  shown  that  if  used  together  with  a  TTM  model  for  convection 
controlled  solid-liquid  phase-change  problems,  this  will  result  in  a 
serious  inconsistency.  This  is  because  TTM  uses  a  mushy  region  to 
guarantee  that  the  temperature  is  continuous  while  in  SOM, 
velocities  are  discontinuous  at  the  solid-liquid  interface.  To  avoid 
this  discontinuity,  a  variant  of  this  method  uses  a  ramped  switch- 
off  of  velocities  by  introducing  a  mushy  zone  between  solid  and 
liquid  phases. 

Instead  of  modifying  coefficients  of  momentum  and  velocity- 
correction  equations,  an  alternate  approach  modifies  the  source 
term  (STM)  in  the  energy  equation  at  the  interface.  As  SOM,  this 
method  also  presents  a  discontinuity  at  the  interface  and, 
therefore,  suffers  from  the  same  problem.  Therefore,  a  ramped 
version  of  this  approach  should  also  be  used  (here  RSTM).  It  should 
be  noted  that  Darcy  STM  developed  for  phase-change  simulations 
in  the  context  of  enthalpy  method  is  a  type  of  ramped  STM  [20,30]. 
Finally,  a  third  approach  is  to  work  with  a  variable  viscosity  of  the 
medium  (WM)  as  proposed  by  Gartling  [31].  In  this  method, 
viscosity  is  set  to  a  very  high  value  in  the  solid  phase,  which 
effectively  stops  all  motions. 

Ma  and  Zhang  [18]  did  a  comparison  between  each  of  the  TTM 
methods  and  the  experimental  results  of  Okada  [32].  They 
concluded  that  ramped  SOM  and  ramped  STM  performed  better 
at  a  slightly  higher  computational  cost.  In  addition,  they  pointed 
out  that  if  the  time  step  is  too  small,  the  convection  is  not  well 
modeled  and  that  the  solutions  of  the  model  diverge. 


h=  f  PkckdT.  (4) 

J  Tm 

In  the  case  of  isothermal  phase  change,  the  liquid  fraction/is  given 
by 

(0  T  <Tm  solid, 

/  =  <  ]0, 1  [  T=Tm  mushy,  (5) 

[  1  T>Tm  liquid. 

Using  Eqs.  (3)  and  (4),  one  can  write  an  alternate  form  for 
dimensional  heat  transfer  in  the  PCM  in  two  dimensions: 


dh 

9 1 

3-  9  1 

dt~ 

dx  ' 

K  dx) 

v  9 y) 

-  M 


df 

dt 


(6) 


where  a  is  the  thermal  diffusivity. 

Very  early,  two  approaches,  finite  difference  and  finite-element 
techniques,  were  used  to  solve  the  phase-change  problems 
numerically.  Then,  finite  volume  and  control  volume  finite- 
element  methods  were  also  employed  [26].  Since  the  introduction 
of  finite-element  and  finite-volume  methods  in  the  1970s,  they  are 
mostly  preferred  over  finite-difference  methods. 

Since  temperature  can  oscillate  with  time  in  some  circum¬ 
stances,  some  authors  have  relayed  on  a  temperature  based  heat 


3.2.  Adaptive  mesh 

To  keep  a  reasonable  representation  of  the  physical  processes 
that  occur  at  the  melt  front,  model  grid  density  must  be  high 
enough  to  smoothly  cover  the  solid-liquid  interface.  However, 
this  high  density  is  not  needed  elsewhere  in  the  numerical 
domain.  Therefore,  it  is  natural  to  adapt  the  grid  density  to  the 
local  physical  conditions  to  improve  the  computational  efficiency. 
There  are  two  main  approaches  used  to  do  this.  One  uses  a  local 
mesh  refinement  method,  i.e.  h-method  [33-36].  In  this  case,  the 
model  starts  with  a  uniform  grid  and,  at  each  iteration,  grid  points 
are  added  or  removed  to  match  the  required  accuracy.  This  is  the 
method  used  in  most  commercial  codes.  The  main  problem  with 
this  method  consists  in  maintaining  data  structures  since  the 
topology  and  the  number  of  grid  points  is  changing  between  each 
time  step. 

To  avoid  this  problem,  the  r-method  (r  stand  for  relocation), 
also  known  as  the  moving  mesh  method,  starts  with  a  uniform  mesh 
and  then  moves  the  mesh  points,  keeping  the  mesh  topology  and 
number  of  mesh  points  fixed  as  the  solution  evolves.  Grid 
deformation  is  usually  done  by  tracking  a  rapid  variation  of  either 
the  solution  or  one  of  its  higher  order  derivatives.  This  method  has 
been  used  with  phase-change  problems  as  in  [37-40]. 
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Lacroix  and  Voller  [41  ]  performed  a  study  comparing  methods 
of  simulation  of  a  phase-change  model  in  a  rectangular  cavity.  They 
concluded  that  the  fixed  grid  must  be  finer  for  a  material  with  a 
unique  melting  temperature  while  the  limiting  factor  in  moving 
mesh  is  the  need  to  use  a  coordinate  generator  at  each  time 
increment.  Comparison  between  fixed  and  moving  grid  has  also  be 
done  by  Viswanath  and  Jaluria  [42]  and  by  Bertrand  et  al.  [43].  In 
this  last  paper,  it  was  found  that  front-tracking  methods  are  better 
adapted  to  the  problem  than  the  fixed-grid  procedures.  The  front¬ 
tracking  methods  would  however  fail  to  simulate  situations  where 
the  transition  from  the  liquid  to  the  solid  phase  is  not  a 
macroscopic  surface,  and  enthalpy  methods  are  to  be  used  in 
most  solidification  problems  where  a  solid-liquid  interfacial 
region  is  present  between  both  phases. 

3.3.  Fist  law  and  second  law  models 


Some  models  are  designed  to  measure  the  first  law  efficiency 
while  others  are  for  the  second  law.  First  law  models  have  some 
shortcomings  because  they  do  not  consider  the  effect  of  time 
duration  through  which  the  heat  is  stored  or  retrieved,  the 
temperature  at  which  the  heat  is  supplied,  and  the  temperature 
of  the  surroundings.  Second  law  models  address  those  issues  which 
lead  to  optimal  designs  and  operations.  Nevertheless,  second  law 
models  are  intended  to  complement  not  replace  the  first  law  models. 

The  first  law  efficiency  of  a  thermal  energy  storage  system  is 
defined  by 


where  Qr  is  the  total  thermal  energy  extracted  from  the  storage 
material  during  the  heat  recovery  process  and  (L  is  the  total  heat 
stored  by  the  phase-change  material  during  the  heat  storage 
process. 

The  second  law  efficiency  of  the  thermal  energy  storage  system 
is  defined  by  the  following  equation: 


iff 


(8) 


where  <pr  is  the  availability  recovered  from  the  storage  system  during 
the  heat  recovery  process  and  <ps  is  the  availability  added  to  the 
system  during  the  heat  storage  process.  This  shows  that  the  first  law 
efficiency  gauges  how  energy  is  utilized,  whereas  the  second  law 
efficiency  indicates  how  well  the  availability  of  energy  is  used. 

An  important  feature  in  the  second  law  analysis  of  the  thermal 
energy  storage  system  is  the  calculation  of  the  entropy  generation 
number  Ns.  This  number  is  defined  as  the  degraded  exergy  divided 
by  the  total  input  exergy  to  the  cycle: 


Ns  = 


Wd 

Ws  +  Wr ' 


(9) 


where  Wd  is  the  total  exergy  destroyed  throughout  the  cycle,  while 
Ws  and  Wr,  are  the  availability  of  the  working  fluid  entering  the 
system  during  the  heat  storage  and  removal  processes,  respective¬ 
ly.  The  second  law  efficiency  can  be  related  to  the  number  of 
entropy  generation  units  by  the  following  equation: 


V/  =  1  -  N,.  (10) 

In  this  situation,  the  system  achieves  its  maximum  second  law 
efficiency  when  Ns  is  zero. 


4.  Model  collection 


In  this  section,  we  review  the  scientific  literature  on  numerical 
models  of  PCM  used  in  thermal  storage  applications.  This  review  is 


based  on  compilations  by  Verma  et  al.  [8]  and  Regin  et  al.  [4] 
combined  with  additional  articles  found  trough  a  recent  literature 
review.  These  papers  have  been  first  classified  by  their  geometrical 
configuration  for  more  convenience  and  then  by  applications. 

4.1.  Rectangular  geometry 

In  a  pioneer  work,  Shamsundar  and  Sparrow  [44,45]  applied  the 
finite-difference  method  to  the  resolution  of  the  enthalpy  equation 
in  the  solidification  of  a  flat  plate.  Shamsundar  also  used  this 
method  to  the  case  of  a  square  geometry  [46-48]. 

Hamdam  and  Elwerr  [49]  presented  a  simple  model  where  all 
sides  of  a  rectangular  enclosure  were  well  insulated  except  one 
vertical  side  where  heat  was  applied.  In  this  model,  the  rate  of 
melting  depends  essentially  on  the  properties  of  the  material,  such 
as  thermal  diffusivity,  viscosity,  conductivity,  latent  heat  of  fusion 
and  specific  heat.  The  model  follows  a  two-dimensional  melting 
process  of  a  solid  PCM  considering  convection  as  the  dominant 
mode  within  the  melt  region,  except  within  the  layer  very  close  to 
the  solid  boundary  in  which  only  conduction  is  assumed  to  take 
place.  Comparison  between  this  model  and  results  of  Bernard  et  al. 
[50]  shows  a  good  agreement  except  that  position  and  inclination 
of  the  melt  front  deviate  as  time  progresses. 

Lacroix  et  al.  [51-54]  solved  the  problem  of  fusion  in  a 
rectangular  cavity  including  the  natural  convection  effects  using  a 
methodology  similar  to  the  moving  mesh  technique.  Results  of 
these  simulations  indicate  that  the  melted  fraction  from  close- 
contact  melting  at  the  bottom  of  the  cavity  is  larger,  by  an  order  of 
magnitude,  than  that  from  the  conduction  dominated  melting  at 
the  top.  The  melting  process  is  essentially  governed  by  the 
magnitude  of  the  Stefan  number.  It  is  also  strongly  influenced  by 
lateral  dimensions  of  the  cavity.  The  importance  of  the  convection 
also  implies  that  system  efficiency  is  increased  when  the  heat  is 
concentrated  on  few  spots  instead  of  being  spread  over  the  whole 
surface.  Later  Brousseau  and  Lacroix  [55,56]  studied  performances 
of  a  multi-layer  PCM  storage  unit.  Their  model  is  based  on  the 
conservation  equation  of  energy  for  the  PCM  and  for  the  HTF  (heat- 
transfer  fluid).  This  model  shows  that  cyclic  melting  and 
solidification  yields  complex  isotherm  distributions  within  the 
PCM  with  multiple  phase  fronts. 

Works  by  Costa  et  al.  [57]  studied  numerically  a  two- 
dimensional  rectangular  area,  using  energy  equations  in  solid 
and  liquid  phases,  continuity,  momentum,  and  Stefan’s  equation  at 
the  boundary.  They  applied  the  enthalpy-porosity  approach  used 
by  Voller  et  al.  [30,250]  to  the  governing  equations.  Three  PCMs 
were  analyzed:  paraffin  (n-octadecanol)  and  metals  (gallium  and 
tin).  They  compared  their  results  with  those  obtained  in  the 
literature  [245-249].  In  the  case  of  octadecanol,  there  is  poor 
agreement  with  experimental  results  in  the  upper  zone  where  melt 
liquid  from  the  sides  fills  the  upper  empty  cavity  and  accelerates 
melting  in  this  area.  Variation  of  viscosity  with  temperature, 
heating  losses  in  the  wall  or  different  initial  supercooling  are 
proposed  explanations  for  this  discrepancy  [57,58]. 

Since  most  PCMs  have  low  thermal  conductivities,  then  heat- 
transfer  rates  limit  their  applications.  To  improve  their  perfor¬ 
mances  they  are  often  used  in  a  thin  plate  configuration  like  a  heat 
exchanger  [59-64].  Costa  et  al.  [65]  also  studied  a  thermal  storage 
system  consisting  of  seven  thin  rectangular  containers  of  PCM.  The 
shape  of  the  container  is  used  to  enhance  the  rate  of  heat  transfer. 
The  performance  of  this  system  has  been  analyzed  using  an 
enthalpy  formulation  and  a  fully  implicit  finite-difference  method 
[25,66].  Vakilaltojjar  and  Saman  [67]  proposed  to  improve  on  this 
design  by  using  PCM  with  different  melting  point  (CaCl2-6H20  and 
potassium  fluoride  tetrahydrate  (I<F-4H20). 

Dolado  et  al.  [68]  developed  a  model  for  thermal  energy  storage 
unit  for  air  conditioning  application  using  thin  PCM  slab.  In 
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parallel,  they  have  also  characterized  the  PCM  in  laboratory  as  well 
as  experimentally  validated  their  model.  Comparing  various 
numerical  approaches,  they  showed  that  a  simple  ID  implicit 
finite  differences  model  seems  to  be  an  appropriate  model  to 
simulate  both,  single  PCM  slabs  and  full  storage  systems. 
Comparing  the  advantages  of  various  numerical  approaches,  they 
concluded  that  semi-analytical  models  give  a  fast  first  approach, 
successful  to  optimize  the  design,  while  2D  model  should  be  used 
to  study  thicker  plates  and  that  the  complete  computational  fluid 
dynamics  model  should  be  employed  for  more  rigorous  analyses. 

Zukowski  [69]  has  analyzed  the  heat  and  mass  transfer  in  a 
ventilation  duct  filled  with  encapsulated  paraffin  wax  RII-56.  This 
author  proposed  a  new  approach  for  approximating  the  specific  heat 
of  the  PCM  as  a  function  of  its  temperature  for  the  whole  range  of 
operating  conditions,  by  using  an  interpolating  cubic  spline  function 
to  estimate  specific  heat  of  PCM  at  each  temperature.  His  analysis 
showed  that  the  average  charge  time  (melting)  of  the  tested  units 
was  over  twice  higher  than  the  time  of  the  discharge  cycle  (freezing). 

Silva  et  al.  [70]  have  modeled  an  enclosure  filled  with  wax  in  a 
rectangular  geometry  oriented  vertically.  Their  results  show  that  a 
simplified  numerical  model  can  be  used  to  predict,  with  reasonable 
accuracy,  the  performance  of  this  kind  of  heat  exchange  unit. 
Vynnycky  and  Kimura  [71]  produced  an  analytical  and  numerical 
study  of  coupled  transient  natural  convection  and  solidification  in 
a  rectangular  enclosure.  Their  work,  based  on  a  non-dimensional 
analysis,  indicates  that  some  asymptotic  simplification  is  possible 
for  materials  commonly  used  in  literature  (water,  gallium,  lauric 
acid).  This  suggests  that  problems  might  be  simplified,  when 
Ra  »  1  and  St  <s  1,  giving  a  conventional  boundary  value  problem 
for  the  liquid  phase  and  pointwise-in-space  first-order  ordinary 
differential  equations  for  the  evolution  in  time  of  the  solidification 


front.  The  method  was  tested  against  full  2D  finite-element-based 
transient  numerical  simulations  of  solidification.  The  asymptotic 
analysis  decouple  the  fluid  flow  and  heat-transfer  problem  in  the 
liquid  from  the  heat-transfer  problem  in  the  solid,  and  is  able  to 
describe  the  quantitative  features  of  the  numerical  solutions  very 
well  for  about  90%  of  the  height  of  the  enclosures.  However,  in  the 
final  10%,  the  analytical  solution  is  not  uniformly  valid  for  all  time, 
and  tends  to  overestimate  the  final  thickness  of  the  solid  layer. 

Zivkovic  and  Fujii  [72]  have  created  a  simple  computational 
model  for  isothermal  phase  change.  The  model  was  based  on  an 
enthalpy  formulation  with  equations  constructed  in  such  a  form 
that  the  only  unknown  variable  is  the  PCMs  temperature.  The 
theoretical  model  was  verified  with  a  test  problem  and  an 
experiment  performed  by  authors.  Experimental  and  computa¬ 
tional  data  demonstrated  that  heat  conduction  within  the  PCM  in 
the  direction  of  the  fluid  flow,  the  thermal  resistance  of  the 
container’s  wall,  and  the  effects  of  natural  convection  within  the 
melt  could  be  ignored  for  the  conditions  investigated  in  this  study. 
This  model  also  shows  that  the  rectangular  container  requires 
nearly  half  of  the  melting  time  as  for  the  cylindrical  container  with 
the  same  volume  and  heat-transfer  area.  Using  the  same 
formalism,  Najjar  and  Hasan  [73]  built  a  model  of  a  greenhouse 
using  a  PCM  as  a  heat  storage  unit.  They  validate  their  model  by 
comparison  with  the  experimental  results  of  [72], 

Other  studies  have  also  explored  the  behavior  of  PCM  in 
rectangular  geometries  [74,75].  They  are  not  explicitly  reviewed 
herein  but  provided  for  the  interested  reader.  Table  2  presents  a 
summary  of  numerical  methods  used  in  the  context  of  the 
Cartesian  coordinates  system.  The  symbols  FG  stand  for  fixed  grid, 
FD  for  finite  difference,  FE  for  finite  element,  A  for  analytic,  CFD  for 
computed  fluid  dynamics  and  MM  for  moving  mesh  respectively. 


Table  2 

Models  for  rectangular  geometry. 


Model 

Material 

Numerical 

Comment 

Validation 

formulation 

Shamsundar  and 

FD 

2D 

Surface  integrated  heat-transfer  rates;  boundary  temperatures 

Sparrow  [44,45] 

solidified  fraction  and  interface  position  all  as  function  of  the  time 

Hamdam  and 

n-Octadecane 

A 

ID 

Propagation  and  inclination  of  the  interface  and  energy  storage  rate 

Experimental  [50] 

Elwerr  [49] 

predicted 

Numerical  [214,215] 

Lacroix  et  al.  [51-54] 

n-Octadecane  [51] 

FG 

2D 

The  melting  process  is  chiefly  governed  by  the  magnitude  of  the 

Experimental  [53,246] 

Stefan  number 

Numerical  [16] 

Brousseau  and 

FD 

2D 

Multi-layer  PCM  storage 

Numerical  [251] 

Lacroix  [55,56] 

FG 

Parametric  model  formulated  from  numerical  study 

Costa  et  al.  [57] 

n-Octadecanol 

FD 

2D 

For  n-octadecanol,  there  is  poor  agreement  with  experimental 

Experimental 

gallium,  tin 

results.  Variation  of  viscosity  with  temperature,  as  heating 
losses  in  the  wall  or  different  initial 

supercooling  are  the  proposed  explanation  for  this  discrepancy 

[245-249] 

Costa  et  al.  [65] 

Multiple  PCM 

FD 

ID 

For  1 D  model,  calculations  have  been  made  for  the  melt 

Numerical 

2D 

2D 

fraction  and  energy  stored  for  conduction  and  convection, 

while  for  2D  model,  calculations  have  been  made  for  conduction  only. 

[216,217] 

Vakilaltojjar  and 

CaCl2-6H20 

A 

Thin  flat  containers  and  air  is  passed  through  gaps  between  them 

On  going,  not 

Saman  [67] 

KF-4H20 

FE 

published 

Dolado  et  al.  [68] 

A 

ID 

Encapsulated  PCM  slab 

Experimental  [218] 

FD 

2D 

CFD 

Zukowski  [69] 

Paraffin  wax 

FD 

3D 

Interpolating  cubic  spline  function  method  is  used  for  determining  an 

Yes 

RII-56 

effective  specific  heat 

Silva  et  al.  [70] 

Paraffin  wax 

MM 

ID 

Charge  and  discharge  of  a  PCM  encapsulated  slab 

Yes 

Vynnycky  and 

Lauric  acid 

MMLE* 

2D 

Laplace-Euler:  new  numerical  method* 

Experimental  [116] 

Kimura  [71] 

Applicable  to  water  and  gallium 

Numerical  [219J 

Zivkovic  and 

CaCl2-6H20 

FD 

ID 

Isothermal  phase  change  of  encapsulated  PCM 

Yes 

Fujii  [72] 

Numerical  [25] 
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4.2.  Spherical  geometry 

Spherical  geometry  represents  an  interesting  case  for  heat 
storage  application,  since  spheres  are  often  used  in  packed  beds  as 
those  described  in  Section  4.4.  Due  to  the  complexity  of  such 
systems,  it  is  often  more  efficient  to  first  model  the  behavior  of  an 
individual  sphere,  to  then  describe  it  with  a  simple  parametric 
model  to  be  used  in  the  packed  bed  modeling. 

Roy  and  Sengupta  [76]  have  examined  the  melting  process  with 
the  solid  phase  initially  uniformly  supercooled.  The  presence  of 
supercooling,  forced  them  to  modify  the  heat-transfer  equation  to 
include  the  effects  of  a  temperature  gradient  in  the  solid  core.  As  a 
result,  a  closed-form  solution  could  not  be  obtained.  At  every  time 
step,  the  unsteady  conduction  equation  has  been  solved  numeri¬ 
cally  using  a  toroidal  coordinate  system,  which  has  been  suitably 
transformed  to  immobilize  the  moving  boundary  and  to  transform 
the  infinite  domain  into  a  finite  one.  In  a  subsequent  paper  [77], 
they  have  examined  the  impact  of  convection  on  the  melting 
process.  They  demonstrated  that  the  fluid  flow  in  the  upper  liquid 
region  is  essentially  quasi-steady,  since  the  liquid  velocity  in  both 
the  film  and  the  upper  zone  are  much  greater  than  the  rate  of 
movement  of  the  interface.  They  also  demonstrated  the  impor¬ 
tance  of  the  ratio  of  the  buoyancy  force  due  to  density  variation  of 
the  liquid  with  temperature,  and  the  density  difference  between 
the  two  phases,  which  reduces  the  melt  rate  and  places  an  upper 
bound  for  the  range  over  which  the  film  solution  is  valid. 

Ettouney  et  al.  [78]  experimentally  evaluated  the  heat  transfer 
in  paraffin  inside  spherical  shell.  Their  calculation  shows  that  the 
Nusselt  number  in  the  phase-change  material  during  melting  is 
one  order  of  magnitude  higher  than  during  solidification  and  is 
strongly  dependent  on  the  sphere  diameter.  In  consequence,  the 
melting  was  about  a  factor  two  slower  than  the  freezing.  This 
counter-intuitive  behavior  is  caused  by  the  convection,  which 
brings  a  large  fraction  of  the  heat  to  the  upper  part  of  the  sphere, 
where  it  heats  the  melt  fraction  instead  of  the  solid  phase.  In  the 
freezing  process,  conduction  dominates  and  the  freezing  occurs 
around  the  entire  surface.  As  the  melt  volume  decreases  and  the 
role  of  natural  convection  diminishes  rapidly.  This  result  contrasts 
with  the  work  of  Barba  and  Spriga  [79]  on  the  behavior  of 
encapsulated  salt  hydrates,  used  as  latent  energy  storage  in  a  heat- 
transfer  system  of  a  domestic  hot  water  tank,  who  found  that  the 
shortest  time  for  complete  solidification  is  provided  by  small 
spherical  capsules,  with  high  Jacob  numbers  and  thermal 
conductivities.  On  the  impact  of  convection,  Khodadadi  and  Zhang 
[80]  also  noted  that  it  created  an  asymmetry:  melting  in  the  upper 
region  of  the  sphere  much  faster  than  in  the  lower  region  due  to  the 
enhancement.  Their  computational  findings  were  verified  through 
qualitative  constrained  melting  experiments  using  a  high-Prandtl 
number  wax  as  the  phase-change  material. 

A  simple  analytical  solution  was  derived  by  Fomin  and  Saitoh  [81  ] 
to  describe  close-contact  melting  within  a  spherical  capsule.  Its 
accuracy  is  estimated  to  be  10-15%  by  comparison  to  numerical 
analysis.  This  tool  is  likely  to  be  useful  when  evaluating  thermal 
energy  storage  systems  which  contain  thousands  of  capsules.  A 
similar  analysis  has  been  done  for  elliptical  capsule.  It  was  found  that 
elliptical  capsules  show  higher  rate  of  melting  than  circular  ones. 
Prolate  capsules  provide  more  effective  melting  than  oblate  ones, 
even  though  they  have  the  same  aspect  ratios  and  vertical  cross- 
sectional  areas  [82,83].  In  addition,  Adref  and  Eames  [84,85]  have 
been  able  to  derive  semi-empirical  equations  that  allow  the  mass  of 
ice  within  a  sphere  to  be  predicted  during  the  freezing  or  melting 
processes.  Complementary  studies  on  quasi-static  models  and  on 
parametric  models  have  been  done  by  others  research  teams  [86-90]. 

Ismail  and  Henriquez  [88]  created  a  parametric  model  for  the 
study  of  the  solidification  of  a  PCM  in  a  spherical  shell.  This  model 
was  based  on  pure  conduction  in  the  PCM  subject  to  boundary 


conditions  of  constant  temperature  or  convection  heat  transfer  on 
the  external  surface  of  the  spherical  shell.  The  model  was  validated 
by  comparison  with  available  similar  models  [91-94]  and  the 
agreement  was  found  to  be  satisfactory.  Like  other  groups,  they 
found  that  the  increase  of  the  size  of  the  sphere  leads  to  increasing 
the  time  for  freezing.  The  same  team  performed  a  specific 
numerical  and  experimental  study  on  water  solidification  [89]. 
The  governing  equations  of  the  problem  and  associated  boundary 
conditions  were  formulated  and  solved  using  a  finite-difference 
approach  and  a  moving  grid  scheme.  A  shortcoming  of  the  model 
was  its  incapacity  of  adequately  handling  supercooling.  Neverthe¬ 
less,  they  conclude  that  shell  material  affects  the  solidification 
time.  Freezing  was  faster  for  metallic  material  (Cu,  Al)  followed  by 
polyethylene,  acrylic  and  PVC.  However,  the  time  difference  was 
rather  small,  which  still  makes  non  metallic  capsules  attractive. 

Veerappan  et  al.  [95]  did  a  an  investigation  the  phase-change 
behavior  of  65  mol%  capric  acid  and  35  mol%  lauric  acid,  calcium 
chloride  hexahydrate,  n-octadecane,  n-hexadecane,  and  n-eico- 
sane  inside  spherical  enclosures.  They  compared  their  numerical 
model  results  with  the  experimental  work  of  Eames  and  Adref  [96] 
and  they  concluded  with  a  good  agreement  between  them, 
notwithstanding  the  fact  that  Eames  and  Adref  [96]  worked  with 
ice  which  floats  unlike  solid  phase  of  modeled  PCMs.  Based  on  their 
modeling,  they  found  that  spheres  with  larger  diameters  take  more 
time  to  completely  solidify  than  smaller  diameters  ones,  since  heat 
has  to  travel  a  smaller  distance  and  hence  rate  of  solidification  is 
higher  for  smaller  capsules.  The  same  is  true  for  melting.  However, 
both  freezing  and  melting  took  approximately  the  same  time, 
which  contrasts  with  result  of  [69,78]. 

Regin  et  al.  [97]  carried  out  both  an  experimental  and  numerical 
analysis  of  a  spherical  capsule  placed  in  a  convective  environment 
filled  with  wax  as  a  PCM.  Model  results  were  compared  to 
experimental  results  and  were  found  to  be  in  good  agreement. 
From  the  model  results,  they  found  that  the  shorter  time  for 
complete  melting  and  higher  time  averaged  heat  flux  occurs  in  the 
capsules  of  smaller  radii  and  for  higher  Stefan  number. 

Lin  and  Jiang  [86]  developed  an  improved  quasi-steady  analysis 
model  for  freezing  in  a  slab,  a  cylinder,  and  a  sphere.  Their 
improvement  was  the  introduction  of  an  additional  term  to  the 
temperature  profile  to  simulate  the  transient  effect  on  the 
temperature  distribution  in  the  solid  phase.  This  additional  term 
is  based  on  the  ratio  of  the  heat  flux  at  the  phase  boundary  to  that 
at  the  cooling  surface,  and  physically,  represents  the  thermal 
capacity  effect  in  the  frozen  region.  The  maximum  relative  error  of 
the  moving  phase  front  location  obtained  from  the  improved 
quasi-steady  analysis  of  a  slab  is  about  3%  in  comparison  with  that 
obtained  from  the  exact  solution  of  the  freezing  process  in  a  slab, 
while  the  error  for  the  classical  quasi-steady  analysis  is  20%.  Since 
there  is  no  exact  solution  available  for  the  freezing  process  taking 
place  in  a  cylinder  or  a  sphere,  the  results  obtained  from  the 
improved  quasi-steady  analysis  are  compared  with  results  of 
[92,98-102]  compiled  by  Lunardini  [103].  For  these  geometries, 
the  maximum  relative  errors  of  the  improved  quasi-steady 
analysis  are  less  than  4%,  while  the  maximum  relative  errors  of 
the  quasi-steady  approximation  are  higher  than  42%.  In  addition  to 
this  improvement  in  precision,  the  improved  model  is  also  less 
sensitive  to  the  variation  of  Stephan  numbers. 

Bilirand  liken  [87]  have  investigated  the  solidification  problem 
of  a  phase-change  material  encapsulated  in  a  cylindrical/spherical 
container.  The  governing  non-dimensional  form  of  equations  of  the 
problem  and  boundary  conditions  were  formulated  and  solved 
numerically  by  using  enthalpy  method  with  control  volume 
approach.  The  numerical  code  was  validated  through  comparison 
with  results  of  other  authors  [88,101,104].  In  general,  there  is  a 
good  agreement.  The  small  differences  are  due  to  the  omission  of 
the  thermal  resistance  of  the  container  and  the  time  dependency  of 
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the  convection  coefficient.  Since  convection  increases  with  the  size  of 
the  container,  differences  get  larger  as  the  sphere  diameter  increases. 
Running  the  simulations  for  various  Stefan  and  Biot  numbers  and 
superheat  parameter,  they  produced  a  regression  which  enables 
them  to  determine  the  coefficient  of  a  function  predicting  the 
freezing  time.  The  coefficient  of  correlation  of  this  regression  is  0.996, 
which  makes  it  suitable  for  solving  engineering  problems. 

Assis  et  al.  [90]  produced  a  parametric  investigation  of  melting 
in  spherical  shells.  Their  simulations  incorporated  such  phenome¬ 
na  as  convection  in  the  liquid  phase,  volumetric  expansion  due  to 
melting,  sinking  of  the  solid  in  the  liquid,  and  close-contact 
melting.  Their  model  attempts  to  solve  complete  transient 
conservation  equations  simultaneously  for  solid  PCM,  liquid 
PCM,  and  air,  while  allowing  for  PCM  expansion,  convection  in 
the  fluid  media  (melted  PCM  and  air),  and  solid  phase  motion  in  the 
liquid  in  a  similar  way  as  in  [105].  This  numerical  simulation 
compares  favorably  with  results  of  the  experimental  investigation 
based  on  previous  work  by  the  same  group  [106,107]. 

Khodadadi  and  Zhang  [80]  performed  a  computational  study  of 
the  effects  of  buoyancy-driven  convection  on  constrained  melting  of 
phase-change  materials  within  spherical  containers.  They  found 
that,  early  during  the  melting  process,  conduction  was  dominant.  As 
the  buoyancy-driven  convection  is  strengthened  due  to  the  growth 
of  the  melt  zone,  melting  in  the  top  region  of  the  sphere  is  much 
faster  than  in  the  bottom  region  due  to  the  enhancement  of  the 
conduction  mode  of  heat  transfer.  They  noted  that  the  intensity  of 
natural  convection  is  more  clearly  indicated  by  the  Rayleigh  number 
than  the  Stefan  number,  as  the  flow  pattern  and  overall  convective 
effects  change  markedly  with  the  Rayleigh  number.  These  numerical 
findings  were  backed  by  an  experimental  visualization  of  the 
melting  of  a  commercial  grade  wax  inside  spherical  bulbs. 

Tan  et  al.  [108,109]  carried  out  a  similar  work  on  constrained 
melting  of  phase-change  materials  to  impend  the  fall  of  solid  phase 
at  the  bottom  of  the  container.  Paraffin  wax  (n-octadecane)  was 
constrained,  through  the  use  of  thermocouples,  during  melting 
inside  a  transparent  glass  sphere.  Their  experimental  setup  provided 
detailed  temperature  data  that  were  collected  during  the  course  of 
the  melting  process  along  the  vertical  diameter  of  the  sphere.  They 
pointed  out  that  experimental  setup  underestimates  the  waviness 


and  the  excessive  melting  of  the  bottom  of  the  PCM  compared  to 
numerical  predictions.  They  concluded  that  this  discrepancy  could 
be  caused  by  the  support  structure  required  to  hold  the  sphere,  since 
it  could  have  inhibited  heat  from  reaching  the  bottom  of  the  sphere. 

Table  3  presents  the  synthesis  for  the  analyses  in  spherical 
geometries. 

4.3.  Cylindrical  geometry 

Anica  Tip  [110,111]  analyzed  the  transient  heat-transfer 
phenomenon  during  technical  grade  paraffin  melting  and  solidifi¬ 
cation  in  a  cylindrical  shell.  The  mathematical  model,  formulated 
to  represent  the  physical  problem,  has  been  proved  suitable  to 
treat  both  melting  and  solidification  processes.  It  can  be  concluded 
that  the  selection  of  the  operating  conditions  and  geometric 
parameters  dimensions  depends  on  the  required  heat-transfer  rate 
and  the  time  in  which  the  energy  has  to  be  stored  or  delivered.  In 
consequence,  these  parameters  must  be  chosen  carefully  to 
optimize  thermal  performances  of  the  storage  unit. 

Gong  and  Mujumdar  [1 12]  developed  a  finite-element  model  for 
an  exchanger  consisting  of  a  tube  surrounded  by  an  external  co-axial 
cylinder  made  up  of  PCMs.  This  model  compares  characteristics  of 
two  operation  modes.  In  mode  1,  hot  and  cold  fluids  are  introduced 
from  the  same  end  of  the  tube  whereas  in  mode  2,  hot  and  cold  fluids 
are  introduced  from  different  ends.  Analyses  show  that  the  energy 
charged  or  discharged  in  a  cycle  using  mode  1  is  5.0%  higher  than 
when  using  mode  2.  The  charge/discharge  rate  is  also  faster  when 
using  mode  1  because  the  temperature  difference  between  the  fluid 
and  the  PCM  is  higher  in  the  fluid  inlet  than  in  the  outlet.  The  larger 
the  temperature  difference  is  the  more  deeply  the  phase-change 
interface  penetrates  into  the  PCM  and  the  more  heat  is  stocked.  In 
the  discharge  process,  symmetrical  phenomena  occur. 

El-Dessouky  and  Al-Juwayhel  [113]  presented  a  technique  to 
predict  the  effect  of  different  designs  and  operating  parameters  on 
the  entropy  generation  number  or  the  second  law  effectiveness  of  a 
latent  heat  thermal  energy  storage  system  (LHTES)  where  the  PCM 
is  contained  in  the  annulus  of  a  double-pipe  heat  exchanger.  One  of 
the  conclusions  of  this  analysis  is  that  second  law  effectiveness 
increases  with  the  increase  of  Reynolds’s  number,  specific  heat- 


Table  3 

Models  for  spherical  geometry. 


Model 

Material 

Numerical 

Comment 

Validation 

formulation 

Roy  and  Sengupta  [76] 

MM 

2D 

Treatment  of  convection 

Numerical  [220] 

FD 

Two  zones  model 

Barba  and  Spriga  [79] 

37.5%  NH4N03  +  62.5% 

A 

ID 

Transient  position  of  interface,  temperature  distribution, 

No 

Mg(N03)-6H20 

melting  fraction,  energy  released,  and  duration  of 
complete  solidification 

Fomin  and  Saitoh  [81] 

n-Octadecane 

A 

2D 

Contact  melting  on  unfixed  solid  phase 

Experimental  [221] 
Numerical  [222,223] 

Adref  and  Eames  [84,85] 

Ice 

Capsule  filled  a  80%  with  an  air  cell  on  the  top 

Eames  and  Adref  [96] 

Ismail  and  Henriquez  [89] 

MM  FD 

ID 

Numerical  correlation  relating  the  working  fluid  temperature 
to  the  time  has  been  produced 

Numerical  [91-94] 

Ismail,  Henriquez 

Ice 

MM  FD 

ID 

Analysis  of  the  impact  of  thermal  conductivity  of  the  shell 

Yes 

and  da  Silva  [88] 

material 

Veerappan  et  al.  [95] 

Various  PCM 

A 

2D 

Solidification  and  melting  of  sphere  with  conduction,  natural 
convection,  and  heat  generation 

Experimental  [96] 

Regin  et  al.  [97] 

Paraffin  wax 

FD 

Convective  environment  outside  capsule 

Yes 

Lin  and  Jiang  [86] 

A 

ID 

Quasi-steady  analysis 

Numerical  [92,98-102] 

Bilir  and  liken  [87] 

FG 

ID 

Correlations  which  express  the  dimensionless  total  solidification 
time  of  the  PCM  in  terms  of  Stefan  Number,  Biot  Number  and 
Superheat  Parameter  were  derived 

Numerical  [88,101,104] 

Assis  et  al.  [90] 

Parrafin  wax  RT27 

MM 

2D 

Partly  filled  capsule  with  open  end 

Yes  +  experimental 
[106,107] 

CFD 

Khodadadi  and  Zhang  [80] 

Beewax 

FG 

2D 

Constrained  melting 

Yes 

Tan  et  al.  [108] 

n-Octadecane 

MM  CFD 

2D 

Constrained  and  unconstrained  melting 

Experimental  [109] 
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transfer  area  and  operating  fluid  wall  temperature  in  the  heat 
recovery  process  and  with  the  decrease  of  fluid  inlet  temperature 
in  the  heat  storage  process. 

For  a  solar  water  heater,  Prakash  et  al.  [114]  analyzed  a  storage 
unit  containing  a  layer  of  PCM  at  the  bottom.  Bansal  and  Buddhi 

[115]  did  the  same  thing  for  a  cylindrical  storage  unit  in  the  closed 
loop  with  a  flat  plate  collector.  In  an  effort  to  improve  performances 
of  phase-change  storage  units  for  solar  heating  applications,  Farid 

[116]  has  suggested  using  many  PCM  with  different  melting 
temperatures.  This  idea  has  been  applied  later  by  Farid  and  Kanzawa 
[  1 1 7  ]  and  Farid  et  al.  [  1 1 8  ]  for  a  unit  consisting  of  vertical  tubes  filled 
with  three  types  of  waxes  having  different  melting  temperatures. 
During  heat  charge,  the  air  flows  first  across  the  PCM  with  the 
highest  melting  temperature  to  ensure  continuous  melting  of  most 
of  the  PCM.  The  direction  of  the  air  flow  must  be  reversed  during  heat 
discharge.  Both  the  theoretical  and  experimental  measurements 
showed  an  improvement  in  performances  of  the  phase-change 
storage  units  using  this  type  of  arrangement. 

Jones  et  al.  [119]  performed  well-controlled  and  well-charac¬ 
terized  experimental  measurements  of  melting  of  a  moderate- 
Prandtl-number  material  (n-eicosane)  in  a  cylindrical  enclosure 
heated  from  the  side.  The  melt  front  was  captured  photographi¬ 
cally  and  numerically  digitized.  A  numerical  comparison  exercise 
was  undertaken  using  a  multi-block  finite-volume  method  and  the 
enthalpy  method  for  a  range  of  Stefan  numbers.  Very  good 
agreement  was  obtained  between  the  predictions  and  the 
experiments  for  Stefan  numbers  up  to  0.1807. 

Esen  and  Ayhan  [120]  developed  a  model  to  investigate  the 
performance  of  a  cylindrical  energy  storage  tank.  In  the  tank,  the  PCM 
was  packed  in  cylinders  and  HTF  flowed  parallel  to  it.  The  PCMs 
considered  were  calcium  chloride  hexahydrate  (CaCl2-6H20),  paraf¬ 
fin  wax  (P-116),  Na2S04-10H20  and  paraffin.  The  performance  of 
cylindrical  energy  storage  tank  was  determined  by  computer 
simulations  [121,122]  and  backed  by  experimental  data.  The  results 
show  that  the  stored  energy  becomes  higher  at  a  given  time  as  the 
mass  flow  rate  or  inlet  HTF  temperature  increases  and  for  more  energy 
storage  appropriate  cylinder  wall  materials  and  dimensions  should  be 
selected,  such  as  higher  thermal  conductivity  and  small  radius. 

Jian-you  [123]  produced  a  numerical  and  experimental 
investigation  of  a  thermal  storage  unit  involving  a  triplex 
concentric  tube  with  PCM  filling  in  the  middle  channel,  with 
hot  heat-transfer  fluid  flowing  into  the  outer  channel  during 
charging  process  and  cold  heat-transfer  fluid  flowing  into  the  inner 
channel  during  discharging  process.  A  simple  numerical  method, 
according  to  energy  conservation,  called  temperature  and  thermal 
resistance  iteration  method  has  been  developed  for  the  analysis  of 
PCM  solidification  and  melting  in  the  triplex  concentric  tube. 
Comparison  between  numerical  predictions  and  the  experimental 
data  shows  good  agreement. 


Table  4  presents  a  complete  picture  at  a  glance  of  numerical 
methods  in  cylindrical  geometries. 

4.4.  Packeds  beds 

Since  most  PCMs  possess  a  low  thermal  conductivity,  increasing 
surface/volume  ratio  is  appealing  to  increase  the  heat-transfer 
rate.  This  can  be  done  by  packing  a  volume  with  a  large  number  of 
PCM  capsules.  Saitoh  and  Hirose  [124]  put  in  evidence  this  idea 
both  theoretically  and  experimentally.  However,  this  improve¬ 
ment  comes  at  the  price  of  a  significant  pressure  drop  and  a  higher 
initial  cost. 

Zhang  et  al.  [125]  produced  a  general  model  for  analyzing 
thermal  performances  of  both  melting  and  solidification  process  of 
a  LHTES  composed  of  PCM  capsules.  This  model  can  be  used  to 
analyze  the  instantaneous  temperature  distribution,  the  instanta¬ 
neous  heat-transfer  rate,  and  the  thermal  storage  capacity  of 
LHTES  system.  From  this  modeling  effort,  it  was  found  that  the 
thermal  performance  of  a  given  system  can  be  determined  by  the 
effective  filling  factor  of  a  PCM  capsule.  Since  formulation  of  the 
effective  filling  factor  for  given  geometric  and  flow  conditions  are 
available  in  the  literature  [1,126-128],  this  largely  reduces 
complexity  of  calculation.  This  model  was  validated  with  the 
results  in  literature  for  a  LHTES  composed  of  PCM  spheres.  The 
model  is  quite  general  so  it  can  be  used  to  optimize  and  to  simulate 
the  thermal  behavior  of  LHTES. 

Benmansour  et  al.  [  1 29  ]  developed  a  two-dimensional  numerical 
model  that  predicts  the  transient  response  of  a  cylindrical  packed 
bed  randomly  packed  with  spheres  of  uniform  sizes  and  encapsu¬ 
lating  paraffin  wax  with  air  as  a  working  fluid.  Experimental 
measurements  of  temperature  distribution  compared  favorably  with 
the  numerical  results  over  a  broad  range  of  Reynolds  numbers  and 
showed  that  the  time  at  which  the  bed  reaches  the  thermal 
saturation  state  decreases  significantly  with  the  mass  flow  rate  of  air. 

Bedecarrats  et  al.  [130-136]  worked  on  a  model  of  a  cylindrical 
tank  filled  with  encapsulated  phase-change  materials.  This  model 
had  the  particularity  of  examining  the  impact  of  supercooling,  which 
can  seriously  deteriorate  performances  of  a  LHTES  [131,132]. 
Numerical  results  were  validated  with  an  experimental  test  tank. 
From  both  numerical  and  experiment  results,  it  was  found  that  the 
optimum  charging  mode  is  obtained  in  the  case  of  the  vertical 
position,  where  the  motions  due  to  the  natural  convection,  are  in  the 
same  direction  as  the  forced  convection.  Also,  they  demonstrated 
that  an  increase  in  flow  rate  increases  the  charge  and  the  discharge 
rates.  They  also  recommend  to  slightly  oversize  the  tank. 

Cheralathan  et  al.  [137]  did  also  a  numerical  and  experimental 
study  of  the  impact  of  super  cooling  and  porosity.  They  came  to  the 
conclusion  that  lower  inlet  HTF  temperature  reduces  the  super¬ 
cooling  of  the  PCM  and  the  total  charging  time.  However,  the 


Table  4 

Models  for  cylindrical  geometry. 


Model 

Material 

Numerical 

formulation 

Comment 

Validation 

Tip  [110,111] 

Parafin  wax  RT  30 

FG 

2D 

Melting + solidification 

Yes 

Gong  and  Mujumdar  [112] 

80.5%  LiF+19.5%  CaF2 

FG 

2D 

Co-axial  exchanger 

No 

El-Dessouky  and  Al-Juwayhel  [113] 

Parafin  wax,  CaCl2-6H20 

A 

2D 

Second  law 

No 

Prakash  et  al.  [114] 

M 

2D 

Vertical  cylinder  with  PCM  at  the  bottom 

Yes 

Bansal  and  Buddhi  [115] 

MM 

2D 

Closed  loop  with  a  flat  plate  collector 

Yes 

Farid  and  Kanzawa  [117] 

MM 

2D 

Several  MCPs  with  different  melting 
temperatures 

Yes 

Jones  et  al.  [119] 

n-Eicosane 

MM 

2D 

Multi-block  FVM 

Yes 

Esen  and  Ayhan  [120] 

CaCl2-6H20  paraffin,  Na2SO410H2O 
paraffin 

FD  FG 

2D 

Energy  storage  reservoir 

Experimental 

[224] 

Jian-you  [123] 

n-Hexacosane 

ID 

Triplex  concentric  tube 

Temperature  and  thermal  resistance 
iteration  method 

Yes 
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reduction  of  heat-transfer  fluid  temperature  also  reduces  the 
performance  of  the  refrigeration  system.  Therefore,  there  is  a  trade 
to  do  between  inlet  temperature  and  overall  performances.  Low 
porosity  increase  storage  capacity  of  the  system  at  the  expense  of  a 
slower  charging/discharging  rate. 

Arkar  and  Medved  [138]  also  studied  a  cylindrical  container 
containing  spheres  filled  with  paraffin.  Their  model  was  adapted  to 
take  into  account  the  non-uniformity  of  porosity  and  the  fluid 
velocity.  Both  are  the  consequence  of  a  small  tube-to-sphere 
diameter  ratio  of  their  system.  A  comparison  of  the  numerical  and 
experimental  results  confirmed  the  important  role  of  PCMs 
thermal  properties  especially  during  slow  running  processes. 

Wei  et  al.  [139]  explored  the  impact  of  geometry  (sphere, 
cylinder,  plate  and  tube)  of  capsules  on  packed  bed  performances. 
Effects  of  the  capsule  diameter  and  shell  thickness  and  the  void 
fraction  were  also  investigated.  They  discovered  that  the  heat 
release  performance  decreases  in  the  order  from  sphere,  to  cylinder, 
to  plate  and  finally  to  tube.  Working  on  another  aspect  of  system 
optimization,  Seeniraj  and  Lakshmi  Narasimhan  [140]  study 
considered  the  impact  of  heat  leak  through  sidewalls  by  convection. 
Their  results  indicate  a  significant  delay  in  the  solidification  process 
with  a  relatively  lower  charging  rate  for  the  LHTS  unit  with 
convective  heat  leak.  Based  on  this  result,  they  derived  limiting 
values  of  heat  leak  ratio  for  effective  charging  of  these  systems. 

Ismail  and  Henriquez  [141]  have  created  a  simplified  transient 
one-dimensional  model  of  PCM  capsule  stocked  in  a  cylindrical  tank. 
The  convection  present  in  the  liquid  phase  of  the  PCM  is  treated  by 
using  an  effective  heat  conduction  coefficient.  The  solution  of  the 
differential  equations  was  realized  by  the  finite-difference  approxi¬ 
mation  and  a  moving  mesh  inside  the  spherical  capsule.  The 
geometrical  and  operational  parameters  of  the  system  were 
investigated  both  numerically  and  experimentally  and  their  influ¬ 
ence  on  the  charging  and  discharging  times  was  investigated.  Ismail 
and  Stuginsky  [142]  compared  four  basic  groups  of  models  of  packed 
beds:  the  continuous  solid  phase  models,  Schumann’s  model,  the 
single  phase  models  and  the  thermal  diffusion  models.  Numerical 
tests  were  carried  out  to  evaluate  the  time  consumed  by  each  of 
them.  It  was  found  that  these  models  were  essentially  equivalent 
except  for  the  computational  time.  For  example,  model  with  a 
thermal  gradient  inside  the  solid  particles  needed  more  than  twenty 
times  the  computing  resources  required  by  Schumann’s  model. 

Chen  [143,144]  also  developed  a  one-dimensional  porous- 
medium  model.  Predictions  of  this  model  were  compared  with 
the  lumped  capacitance  model  which  assumes  a  uniform  tempera¬ 
ture  in  the  packed  capsules  and  in  the  coolant  flow  and  with 
experimental  data  derived  from  literature.  Results  show  that  the 
one-dimensional  model  has  the  advantage  of  predicting  the 
temperature  distributions  of  PCM  and  coolant.  For  the  study  of 
transient  response  of  PCM  beds,  Goncalves  and  Probert  [145]  have 
introduced  a  new  non-dimensional  number,  the  Gonsalves  number, 
which  mean  integrated  value  can  be  used  to  characterize  the 
thermal  response  of  a  thermal  energy  storing  system. 

Adebiyi  et  al.  [146]  have  developed  a  model  of  a  packed  bed  for 
high-temperature  thermal  energy  storage  that  incorporates  several 
features:  thermal  properties  variations  with  temperature,  provi¬ 
sions  for  arbitrary  initial  conditions  and  time-dependent  varying 
fluid  inlet  temperature,  formulation  for  axial  thermal  dispersion 
effects  in  the  bed,  modeling  for  intraparticle  transient  conduction, 
energy  storage  in  the  fluid  medium,  transient  conduction  in  the 
containment  vessel  wall.  Energy  extraction  can  be  achieved  either 
with  flow  direction  parallel  to  that  in  the  storage  mode  (concurrent) 
or  in  the  opposite  direction  (counter-current).  This  model  also 
allows  the  computation  of  the  first  and  second  law  efficiencies.  Yagi 
and  Akiyama  [147]  have  also  studied  the  heat-transfer  problem  at 
high  temperatures.  Six  different  materials  were  tested  as  PCMs:  two 
inorganic  and  four  metallic  materials.  The  metal  PCMs  were  found  to 


be  excellent  for  heat  storage  because  of  uniform  temperature  in  the 
capsule.  Heat-transfer  simulation  was  also  conducted  for  a  packed 
bed  process  of  spherical  capsules  and  concurrent  flow  for  heat 
storage  and  release  showed  better  result  for  effective  use  of  storage 
heat  than  counter-current  flow. 

Wanatabe  et  al.  [148,149]  have  explored  the  impact  of 
combining  PCMs  with  various  melting  points  on  exergy  efficiency 
and  charging  and  discharging  rates  in  a  latent  heat  storage  system. 
The  heat  storage  module  consisted  of  horizontal  cylindrical 
capsules  filled  with  three  types  of  PCMs.  The  optimum  melting 
point  distribution  of  the  PCMs  has  been  estimated  from  numerical 
simulations  and  also  from  simple  equations.  The  fast  charging  or 
discharging  rate  leads  to  high  exergy  efficiency.  Both  the 
experimental  and  numerical  results  showed  some  improvement 
in  charging  and  discharging  rates  by  use  of  multiple  PCMs. 

MacPhee  and  Dincer  [150]  studied  the  melting  and  freezing  of 
water  in  spherical  capsules  using  ANSYS,  GAMBIT,  and  FLUENT  6.0. 
These  models  were  validated  through  a  rigorous  time  and  grid 
independence  tests  and  their  numerical  results  were  successfully 
compared  to  the  experimental  ones  of  [78].  Energy  efficiency  was 
calculated  and  found  to  be  over  99%,  though  viscous  dissipation  was 
included.  Using  exergy  analysis,  the  exergetic  efficiencies  are 
determined  to  be  about  75%  to  over  92%,  depending  on  the  HTF 
scenario.  The  two  efficiencies  decreased  when  hot  transfer  fluid  flow 
was  increased,  due  mainly  to  increasing  heat  losses  and  exergy 
dissipation.  These  results  are  different  from  those  obtained  by 
Wanatabe  et  al.  [148,149].  Higher  temperature  differences  between 
hot  fluid  and  melting  temperature  of  PCMs  were  found  to  be  most 
optimal  exergetically,  but  least  optimal  energetically  due  to  entropy 
generation  accompanying  heat  transfer.  Higher  temperature 
differences  also  increased  the  rates  of  charging/discharging. 
However,  increasing  the  HTF  flow  rate  did  not  decrease  the  charging 
and  discharging  times  by  any  relatively  significant  amount.  This 
result  is  somewhat  discrepant  with  various  previous  studies  [129- 
136].  The  same  team  also  developed  a  simple  analytical  model  [151  ] 
that  describes  heat  transfer,  fluid  flow,  and  thermodynamic,  which 
brought  essentially  the  same  conclusions  as  the  numerical  model. 

A  summary  for  packed  beds  models  is  provided  in  Table  5. 

4.5.  Finned  geometry 

Another  approach  to  increase  the  heat  exchange  efficiency  in 
PCMs  is  to  use  fins  in  the  heat  exchange  system.  These  somehow 
more  complex  geometries  have  been  extensively  studied. 

Sasaguchi  et  al.  [152,153]  studied  experimentally  and  theoreti¬ 
cally  effects  of  the  configuration  of  a  finned  tube  on  the  heat-transfer 
characteristics  of  a  latent  heat  thermal  energy  storage  system.  They 
concluded  that  performance  of  the  unit  is  almost  the  same  for  any 
configuration  with  the  same  surface  area  even  if  it  is  quite  different. 

Lacroix  [154]  developed  a  theoretical  model,  using  an  enthalpy- 
based  method  coupled  to  a  convective  heat  transfer  of  a  shell-tube 
thermal  storage  unit  with  PCM  on  the  shell  side  and  the  HTF 
circulating  inside  the  tubes.  The  numerical  model  was  validated 
through  comparison  with  experimental  data.  In  addition,  compar¬ 
ison  was  done  between  models  of  bare  and  finned  tubes.  Results 
showed  that  the  annular  fins  are  the  most  effective  for  moderate 
mass  flow  rates  and  small  inlet  temperatures.  Lacroix  and 
Benmadda  [  1 55  ]  also  studied  the  behavior  of  a  vertical  rectangular 
cavity  filled  with  PCM.  They  used  a  fixed-grid  model  which  was 
validated  with  experimental  data.  They  found  that  both  solidifica¬ 
tion  and  melting  rates  were  improved  by  long  fins.  However, 
melting  rates  reach  a  maximum  when  the  number  of  fins  is 
increased  since  they  interfere  with  the  convection. 

Zhang  and  Faghri  [156]  studied  the  heat-transfer  enhancement 
produced  by  externally  finned  tube  in  a  latent  heat  energy  storage 
system.  The  heat  conduction  in  the  tube  wall,  fins  and  the  PCM 
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Model 

Material 

Numerical 

formulation 

Comment 

Validation 

Zhang  et  al.  [125] 

A 

ID 

A  general  model,  great  reference  paper 

Numerical  [225-229] 

Benmansour  et  al.  [129] 

Parafin  wax 

FD 

2D 

Mass  flow  rate  is  determinant 

Yes 

Bedecarrats  et  al.  [130-136] 

PM 

2D 

Supercooling  studied 

Experimental  [131,135,230] 

Cheralathan  et  al.  [137] 

Ice 

ID 

Supercooling  is  reduced  by  lower  inlet  HTF  temperature 

Yes 

Arkar  and  Medved  [138] 

Parafin  wax  RT20 

PM  FD  FG 

2D 

Continuous  solid  phase  model 

Yes 

Wei  et  al.  [139] 

Paraffin  wax  FNP-0090 

FD  FG 

ID 

Geometry  analysis  of  capsules 

Yes 

Seeniraj  and  Lakshmi 

Heat  leaks  through  sidewalls 

Experimental  [144] 

Narasimhan  [140] 

Ismail  and  Henriquez  [141] 

Ice 

FD  MM 

ID 

Effective  heat  conduction  coefficient 

Yes 

Ismail  and  Stuginsky  [142] 

ID 

Parametric  model 

Chen  [143,144] 

Ice 

ID 

Solution  by  Laplace  transform 

Experimental  [144] 

Goncalves  and  Probert  [145] 

MgCl26H20 

Their  own  dimensionless  number! 

Yes 

Adebiyi  et  al.  [146] 

Zr02 

Cu 

KN03-NaN03 

First  and  second  law  efficiency,  high  temp 

Yes 

Yagi  and  Akiyama  [147] 

ID 

Six  materials  with  various  melting  points 

Yes 

NaCl,  Pb 

Al-Si 

Wanatabe  et  al.  [148,149] 

ID 

Second  Law  model 

Yes 

MacPhee  and  Dincer  [150] 

2D 

First  and  second  law  efficiency.  Fluent 

Experimental  [78] 

were  described  by  a  temperature  transforming  model  using  a 
fixed-grid  numerical  method  [27].  This  model  assumes  that  the 
melting  process  occurs  over  a  range  of  phase-change  tempera¬ 
tures,  but  it  also  can  be  successfully  used  to  simulate  the  melting 
process  occurring  at  a  single  temperature  by  taking  a  small  range  of 
phase-change  temperature.  This  model  was  validated  through 
comparison  with  Lacroix  [154]  and  agreement  between  the  two 
was  found  to  be  quite  good.  The  results  show  that  the  local  Nusselt 
number  inside  the  tube  cannot  be  simply  given  by  the  Graetz 
solution.  The  model  also  shows  that  transverse  fins  are  a  more 
efficient  way  to  enhance  the  melting  heat  transfer  if  the  initial 
subcooling  exists  in  the  PCM.  The  height  of  the  fins  has  significant 
effect  on  the  melting  front  on  the  two  sides  of  the  fins  but  has  no 
significant  effect  on  the  melting  front  between  the  transverse  fins. 

Velraj  et  al.  [157]  studied  the  impact  of  internal  longitudinal 
fins  on  a  cylindrical  vertical  tube  filled  with  paraffin.  A  theoretical 
model  that  accounts  for  the  circumferential  heat  flow  through  the 
tube  wall  was  developed  using  enthalpy  formulation  and  was 
employed  in  conjunction  with  the  fully  implicit  finite-difference 
method  to  solve  the  solidification  in  the  convectively  cooled 
vertical  tube.  They  also  generalized  the  enthalpy  formulation  of 
Date  [158]  to  the  case  where  the  material  has  a  range  of  phase- 
change  temperature.  This  model  was  validated  with  experimental 
data.  From  this  analysis,  they  concluded  that  the  solidification  time 
is  approximately  1  /n  time  the  no  fin  case,  where  n  is  the  number  of 
fins.  They  also  pointed  out  that,  in  the  central  part,  the  flow  is 
restricted  between  fins.  Therefore,  when  the  number  of  fins 
exceeds  4,  half  the  fins  can  extend  to  only  half  the  cylinder  radius. 
In  their  subsequent  work  [159],  they  pointed  out  the  necessity  to 
include  the  effect  of  circumferential  heat  flow  through  the  tube 
wall  for  higher  values  of  Biot  numbers  in  order  to  correctly  predict 
the  heat-transfer  behavior.  For  lower  Biot  numbers,  addition  of  fins 
makes  the  surface  heat  flux  more  uniform,  whereas  for  higher  Biot 
numbers  the  addition  of  fins  improves  the  magnitude  of  the 
surface  heat  flux  and  appreciably  reduces  the  solidification  time. 
For  a  given  quantity  of  heat  to  be  extracted,  a  combination  of  lower 
Biot  number  and  higher  Stefan  number  is  recommended  for  the 
uniform  extraction  of  heat.  The  same  team  [160]  also  investigated 
transient  behavior  of  a  finned-tube  latent  heat  thermal  storage 
module  with  thin  circumferential  fins.  Numerical  results  indicate 
an  appreciable  enhancement  in  the  energy  storage  process  with 
the  addition  of  fins  in  the  module. 

Lamberg  and  Siren  [161-164]  developed  an  analytical  model 
for  analysis  of  the  melting  and  solidification  in  a  finned  PCM.  The 


results  of  the  derived  analytical  model  were  compared  to 
numerical  results  of  a  two-dimensional  model  solved  numerically 
by  means  of  the  effective  heat  capacity  method  and  enthalpy 
method.  The  two-dimensional  results  were  compared  to  simplified 
one-dimensional  results  in  order  to  find  out  the  accuracy  of  the 
simplified  analytical  model.  The  results  of  the  experimental  work 
were  then  compared  to  the  numerical  results  calculated  using  the 
effective  heat  capacity  method  and  enthalpy  method,  in  order  to 
find  out  the  accuracy  of  different  numerical  methods.  The  results 
show  that  the  analytical  models  give  a  satisfactory  estimate  of  fin 
temperature  and  of  the  solid-liquid  interface. 

Castell  et  al.  [165]  studied  both  experimentally  and  numerically 
the  impact  of  external  longitudinal  fin  on  PCM  module  placed  in  a 
water  tank  working  with  a  solar  heating  system.  The  geometry  of  the 
PCM  modules  is  the  same  used  in  an  other  experimental  work  done 
by  the  authors  [166,167]  but  with  added  fin  on  the  PCM  module. 
Experimental  work  was  carried  out  to  determine  natural  convection 
heat-transfer  coefficients  for  PCM  modules.  To  compare  the 
behavior  of  the  system  using  finned  PCM  modules  with  another 
one  without  fins,  two  numerical  codes  were  implemented.  The 
results  obtained  from  the  simulation  were  compared  with  the 
experimental  ones  to  validate  the  numerical  codes.  Based  on  those 
experimental  results,  useful  Nusselt  correlations  in  function  of 
Rayleigh  number  were  found  to  evaluate  the  natural  convection 
heat-transfer  coefficient  for  that  specific  geometry. 

Reddy  [168]  modeled  a  PCM-water  solar  system  consisting  of  a 
double  rectangular  cross  section  enclosure  where  the  top  enclosure 
is  filled  with  paraffin  wax  and  the  bottom  is  filled  with  water.  The 
numerical  modeling  has  been  carried  out  using  the  commercial  CFD 
software  FLUENT.  The  geometiy  modeling  and  mesh  were  generated 
in  GAMBIT.  An  enthalpy  porosity  technique  was  used  in  FLUENT  for 
modeling  the  solidification/melting  process.  In  this  technique,  the 
mushy  zone  is  modeled  as  a  pseudo  porous  medium  in  which  the 
porosity  decreases  from  1  to  0  as  the  material  solidifies.  Simulation 
results  indicate  that  an  optimal  number  of  fins  should  be  used  to 
optimize  the  overall  system  performance. 

Gharebaghi  and  Sezai  [169]  investigated  the  enhancement  of 
energy  storage  rate  of  a  thermal  energy  storage  unit  filled  with  a 
PCM  by  inserting  a  fin  array.  Heat  is  transferred  to  the  unit  through 
the  container  walls,  to  which  fins  are  attached;  a  commercial 
paraffin  wax  is  stored  between  the  fins.  A  mathematical  model, 
based  on  the  finite-volume  method  for  a  two-dimensional  domain, 
is  developed  for  solving  the  melting  problem.  A  fixed-grid  method 
is  used  to  simulate  the  melting  of  the  PCM.  The  interface  between 


122 


Y.  Dutil  et  al.  / Renewable  and  Sustainable  Energy  Reviews  15  (201 1)  112-130 


the  solid  and  the  melt  is  modeled  as  a  porous  medium.  Transient 
simulations  have  been  performed  on  non-uniform  grids  for 
different  fin  and  PCM  layer  thicknesses.  For  horizontal  and  vertical 
arrangements  of  modules,  it  was  observed  that  for  high  tempera¬ 
ture  differences,  the  heat-transfer  rate  can  be  increased  as  much  as 
80  times  by  adding  a  fin  array.  The  minimum  enhancement  rate 
was  3  times,  which  was  achieved  with  widely  spaced  fins.  It  was 
also  observed  that  the  Nusselt  number  is  higher  for  vertical 
modules  arrangement  compared  with  that  of  horizontal  modules 
for  all  fin  spacings  and  temperature  differences. 

Ismail  et  al.  [170]  studied  a  thermal  numerical  model  for  the 
solidification  of  PCM  around  a  radially  finned  tube  with  a  constant 
wall  temperature.  Their  model  was  based  on  a  pure  conduction 
formulation  and  the  enthalpy  method.  Numerical  simulations 
were  performed  to  investigate  the  effects  of  the  number  of  fins,  fin 
thickness,  fin  material,  aspect  ratio  of  the  tube  arrangement  and 
the  tube  wall  temperature.  The  study  shows  that  the  geometrical 
aspects  have  strong  influences  on  the  time  for  complete  solidifica¬ 
tion  and  on  the  solidification  rate.  Also,  the  tube  material  as  well  as 
the  tube  wall  temperature  affect  the  time  for  complete  solidifica¬ 
tion.  The  same  team  [171  ],  following  the  same  procedure,  analyzed 
the  behavior  of  the  solidification  of  a  PCM  around  a  vertical  axially 
finned  isothermal  cylinder.  Their  numerical  model  was  validated 
by  experimental  data.  From  this  analysis,  the  authors  suggested 
that  a  metallic  tube  fitted  with  four  or  five  fins  of  constant 
thickness  equal  to  the  tube  wall  thickness  and  of  radial  length 
around  twice  the  tube  diameter  should  be  the  best  compromise 
solution  between  efficiency,  increase  in  the  heat  flow  rate  and  the 
loss  of  available  storage  capacity. 

Kayansayan  and  Ali  Acar  [172]  modeled  the  formation  of  ice 
around  a  finned-tube  heat  exchanger.  Due  to  the  small  tempera¬ 
ture  difference  between  the  tube  outside  surface  and  the  PCM  (less 
than  5  °C),  natural  convection  is  very  weak  and  even  more 
hampered  by  the  presence  of  vertical  fins.  The  validity  of  their 
model  was  established  by  comparison  with  Lacroix  [154]  and  an 
experimental  setup. 

The  thermal  control  of  portable  electronic  devices,  with  their 
intermittent  operation,  is  well  adapted  to  PCM  utilization. 
However,  careful  optimization  needs  to  be  done  to  insure  the 
adequate  sizing  and  consistent  behavior  under  realistic  operation¬ 
al  condition. 

Akhilesh  et  al.  [173]  developed  a  model  of  a  composite  heat  sink 
composed  of  vertical  arrays  of  fins  surrounded  by  phase-change 
material.  It  should  be  noted  that  heat  load  was  coming  from  the  top 
of  the  heat  sink.  This  heat  sink  was  to  be  used  for  cooling  electronic 
package  below  a  set  point  temperature.  This  model  provides  a 
thermal  design  procedure  for  proper  sizing  of  the  heat  sink,  for 
maximizing  energy  storage  and  for  operation  time.  Based  on  a 
scaling  analysis  of  the  governing  two-dimensional  unsteady 
energy  equations,  a  relation  between  the  critical  dimension  for 
the  heat  sink  and  the  amount  of  PCM  used  was  developed. 

Shatikian  et  al.  [105,174,175]  explored  numerically  the  process 
of  melting  of  a  phase-change  material  (PCM)  in  a  heat  storage  unit 
with  internal  fins  open  to  air  at  its  top.  The  phase-change  material, 
paraffin  wax,  was  stored  between  the  fins.  Transient  three  and 
two-dimensional  simulations  were  performed  using  FLUENT  6.0, 
yielding  temperature  evolution  in  the  fins  and  the  PCM.  Using  a 
dimensional  analysis  to  generalize  their  results,  they  have  shown 
that  within  the  same  geometry,  the  Nusselt  number  and  melt 
fraction  depend  on  the  product  of  the  Fourier  and  Stefan  numbers. 
For  relatively  wide  vertical  PCM  layers,  the  Rayleigh  number  must 
also  be  included,  to  take  into  account  the  effects  of  convection  at 
advanced  stages  of  the  melting  process. 

Behunek  and  Fiala  [176]  also  modeled  a  heat  sink  for  electronic 
package  including  PCM  based  on  the  work  of  Nayak  et  al.  [177].  In 
their  configuration,  the  inorganic  crystalline  salt  acting  as  PCM  was 


placed  in  a  chamber  right  over  the  chip.  This  chamber  was  put  in 
contact  with  fin  cooled  by  forced  convection  of  air.  They  first 
created  an  analytical  description  and  a  solution  of  heat  transfer, 
melting  and  freezing  process  in  ID  and  compared  their  solution  to 
the  numerical  solution  obtained,  through  a  FEM  solution  in  ANSYS 
software,  of  a  real  3D  cooler.  In  addition,  results  of  3D  numerical 
solutions  were  verified  experimentally. 

Saha  et  al.  [178]  also  made  an  analysis  of  a  similar  system  and 
geometry  both  experimentally  and  numerically.  Their  experimen¬ 
tal  setup  consisted  of  an  electrical  heater  to  simulate  the  electronic 
chips  with  the  heat  sink  filled  with  eicosane  above.  The  governing 
equations  of  their  numerical  model  follow  a  single  domain 
approach  where  both  fluid  flow  and  heat-transfer  equations  are 
solved  simultaneously.  The  fin,  heater,  and  substrate  materials  are 
assigned  very  high  viscosity  and  melting  points,  which  ensures 
that  they  remain  solid.  In  consequence,  a  common  set  of  governing 
conservation  equations  for  both  the  solid  and  liquid  regions  can  be 
used  [177].  An  optimal  fin  fraction  of  8%  was  found.  This  was 
attributed  to  convection  cells  in  the  molten  region  of  the  PCM 
present  at  a  large  volume  fraction  of  PCM  but  otherwise 
suppressed.  It  was  also  found  that  a  large  number  of  narrow 
cross  section  fins  is  preferable,  since  it  improves  the  thermal 
contact  for  a  same  volume. 

Kandasamy  et  al.  [179]  have  also  studied  a  similar  heat  sink 
experimentally  and  numerically.  Results  show  that  increased 
power  inputs  enhance  the  melting  rate  as  well  as  the  thermal 
performance  of  the  PCM-based  heat  sinks  until  the  PCM  is  fully 
melted.  A  three-dimensional  computational  fluid  dynamics  model 
was  used  to  simulate  the  problem.  To  describe  the  PCM-air  system 
with  a  moving  internal  interface  but  without  inter-penetration  of 
the  two  fluids,  the  volume-of-fluid  model  has  been  used 
successfully  [180]  and  showed  good  agreement  with  experimental 
data.  The  same  team  [181]  using  the  same  approach  explored  the 
effect  of  orientation  on  performances.  They  concluded  that 
orientation  has  a  minimal  effect  in  their  configuration. 

In  an  original  approach,  Wang  et  al.  [  1 82  ]  created  a  constructal 
rule  for  the  design  of  fins  in  PCM  application.  Constructal 
optimization  of  heat  conduction  with  phase-change  mimics  the 
growth  of  plants  in  nature.  The  high  conductivity  material  is 
treated  as  the  plant  root,  and  the  PCM  acts  as  the  soil.  The  product 
of  thermal  conductivity  and  temperature  gradient  integration 
over  time  in  the  whole  melting  time  interval  of  the  PCM  is  taken  as 
the  criterion  to  determine  the  arrangement  position  of  the  high 
conductivity  material.  Its  physical  meaning  is  the  quantity  of  total 
heat  transported  during  the  melting  time  interval.  The  generation 
rule  of  constructal  optimization  of  heat  conduction  with  phase 
change  is  to  grow  the  high  conductivity  material  at  the  position 
with  the  maximum  total  heat  transported  during  the  time 
interval.  The  degeneration  rule  is  to  withdraw  the  high 
conductivity  material  at  the  position  with  minimum  total  heat 
transported  during  the  time  interval.  Comparison  between 
constructural  and  man-made  shapes  confirm  that  constructal 
ones  are  superior  in  melting  time  and  in  the  average  cold  release 
power. 

Table  6  presents  the  synthesis  for  finned  surfaces. 

4.6.  Porous  and  fibrous  materials 

Erk  and  Dudukovic  [183]  modeled  and  experimented  with  a 
PCM  consisting  of  n-octadecane  retained  by  capillary  forces  in  a 
porous  silica  support.  It  optically  changes  from  opaque  to 
semitransparent  when  the  wax  melts,  thereby  allowing  the  melt 
front  within  the  bed  to  be  tracked.  This  property  and  the  outlet 
temperature  were  compared  with  the  model.  Measured  outlet 
temperature  compares  qualitatively  with  model  predictions.  The 
model  quantitatively  predicts  melt  front  movement  in  the  first  60% 
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Models  for  finned  surfaces. 
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Model 

Material 

Numerical 

formulation 

Comment 

Validation 

Sasaguchi  et  al.  [152,153] 

3D 

Performance  independent  of  geometry  only  surface  area 

Yes 

Lacroix  [154] 

n-Octadecane 

FD 

2D 

An  analytic  ID  model  is  also  produce 

Yes 

Lacroix  and  Benmadda  [155] 

FG 

Zhang  and  Faghri  [156] 

FD 

2D 

Validation  of  semi-analytic  result  of  [231] 

Numerical  Lacroix  [154] 

FG 

Kays  and  Crawford  [232] 

Velraj  et  al.  [157,159,160] 

Paraffin  RT  60  [157] 

FD 

2D 

Phase-change  material  kept  inside  a  longitudinal 

Yes  [157] 

LiF-MgF2  [160] 

FG 

internally  finned  vertical  tube 

159]  num.  [155] 

[160]  num.  [233] 

Lamberg  and  Siren  [161-164] 

Paraffin 

FE 

2D 

Analytical  model  validated 

Yes 

Castell  et  al.  [165,234] 

Sodium  acetate 

A 

Horizontal  and  vertical  fins  around  circular  vertical  cylinder 

Yes 

tri  hydrate 

with  graphite 

Reddy  [168] 

Paraffin  wax 

MM 

2D 

Double  rectangular  enclosure  where  the  top  enclosure  is 

Experimental  [252] 

filled  with  paraffin  wax  and  the  bottom  is  filled  with  water 

Gharebaghi  and  Sezai  [169] 

Paraffin  RT27 

FD 

2D 

Slabs  containing  paraffin  and  metal  fins  made  of  aluminum 

Yes 

Ismail  et  al.  [170,171] 

Ice  [105] 

FD 

2D 

Radially  finned  horizontal  tube  [170] 

Yes  [171] 

Paraffin  [174] 

Axially  finned  vertical  tube  [171] 

Kayansayan  and  Acar  [172] 

Ice 

FD 

2D 

Radially  finned  horizontal  tube 

Numerical  Lacroix  [154] 

FG 

Zhang  and  Faghri  [156] 

Akhilesh  et  al.  [173] 

n-Eicosane 

FD 

2D 

Array  of  vertical  fins 

Numerical  Bejan  [235] 

Shatikian  et  al.  [105,174,175] 

Paraffin  wax  RT25 

2D 

Array  of  vertical  fins 

Numerical  [105] 

Behunek  and  Fiala  [176,163] 

CaCl2-6H20 

ID 

3D 

Also  analytical  1 D 

Yes 

of  the  bed.  Discrepancies  between  the  model  and  experiments  are 
linked  to  significant  heat  losses  in  the  small  insulated  bench-scale 
apparatus  used  for  the  experiment.Mesalhy  et  al.  [  1 84,1 85]  studied 
numerically  and  experimentally  the  effect  of  porosity  and  thermal 
properties  of  a  porous  medium  filled  with  PCM.  In  their  model,  the 
governing  partial  differential  equations  describing  the  melting  of 
phase-change  material  inside  porous  matrix  were  obtained  from 
volume  averaging  of  the  main  conservation  equations  of  mass, 
momentum,  and  energy.  Due  to  the  huge  difference  in  thermal 
properties  between  the  phase-change  material  and  the  solid 
matrix,  two  energy  equations  model  was  adopted  to  solve  the 
energy  conservation  laws.  This  model  can  handle  local  thermal 
non-equilibrium  condition  between  the  PCM  and  the  solid  matrix. 
Finite-volume  technique  in  conjunction  with  numerical  grid 
generation  based  on  body  fitted  coordinate  transformation  has 
been  adopted.  From  their  model,  it  was  found  that  the  best 
technique  to  enhance  the  response  of  the  PCM  storage  is  to  use  a 
solid  matrix  with  high  porosity  and  high  thermal  conductivity.  The 
same  team  [186]  investigated  the  properties  of  graphite  foams 
infiltrated  with  phase-change  materials.  The  geometry  studied  was 
two  concentric  cylinders,  the  inner  being  kept  at  a  constant 
temperature  while  the  outer  one  is  filled  with  foam/PCM  matrix. 
The  model  has  been  validated  by  applying  it  on  the  case  described 
by  Khillarkar  et  al.  [187]  and  has  shown  to  be  in  good  agreement. 
Model  results  indicate  that  estimated  value  of  the  average  output 
power  using  carbon  foam  of  porosity  97%  is  about  five  times 
greater  than  that  for  using  pure  PCMs. 

Fukai  et  al.  [188,189]  modeled  and  experimented  the  use  of 
carbon-fiber  combined  with  n-octadecane.  Carbon  fibers  were 
added  to  improve  the  thermal  conductivities  of  phase-change 
materials  packed  around  heat-transfer  tubes.  The  transient 
thermal  responses  improved  as  the  volume  fraction  of  the  fibers 
and  the  brush  diameter  increased.  However,  these  do  not  further 
improve  when  the  diameter  of  the  brush  is  larger  than  the  distance 
between  the  tubes  due  to  thermal  resistance  near  the  tube  walls. 
The  discharge  rate  using  the  brushes  with  a  volume  of  one  percent 
was  about  30%  higher  than  that  using  no  fibers.  In  the  charge 
process,  the  brushes  prevent  natural  convection.  Still,  the  charge 


rate  with  the  brushes  is  10-20%  higher  than  that  with  no  fibers.  A 
heat-transfer  model  describing  anisotropic  heat  flow  in  the 
composite  was  numerically  solved.  The  calculated  transient 
temperatures  agreed  well  with  the  experimental  data.  A  simple 
model  was  also  developed  to  predict  the  heat  exchange  rate 
between  the  composite  and  the  heat-transfer  fluid.  Hamada  et  al. 
[190]  compared  experimentally  and  numerically  the  thermal 
behavior  of  carbon  chip  and  carbon  brush.  The  carbon-fiber  chips 
are  efficient  for  improving  the  heat-transfer  rate  in  PCMs. 
However,  the  thermal  resistance  near  the  heat-transfer  surface 
is  higher  than  that  for  the  carbon  brushes.  As  a  result,  the  overall 
heat-transfer  rate  for  the  fiber  chips  is  lower  than  that  for  the 
carbon  brushes.  Consequently,  the  carbon  brushes  are  superior  to 
the  fiber  chips  for  the  thermal  conductivity  enhancement  under 
these  experimental  conditions. 

Elgafy  and  Lafdi  [191]  studied  experimentally  and  analytically 
the  behavior  of  nanocomposite  carbon  nanofibers  filled  with 
paraffin  wax.  An  analytical  model  was  introduced  based  on  a  one¬ 
dimensional  heat  conduction  approach  to  predict  the  effective 
thermal  conductivity  for  the  new  nanocomposites.  Their  findings 
showed  good  agreement  with  the  experimental  data  of  the  authors 
and  of  Lafdi  and  Matzek  [192]. 

Frusteri  et  al.  [193]  studied  experimentally  and  numerically 
thermal  conductivity  and  charging-discharging  kinetics  of 
inorganic  PCM44  containing  carbon  fibers.  A  numerical  model, 
based  on  the  enthalpy  formulation  and  on  an  implicit  finite- 
difference  method,  has  been  developed  to  predict  the  one¬ 
dimensional  heat-transfer  diffusion  of  a  PCM  containing  carbon 
fibers.  The  results  obtained  by  use  of  this  model  were  satisfactory 
compared  to  experimental  results.  Authors  also  pointed  out  that 
the  values  of  thermal  conductivity  calculated  by  Fukai’s  method 
are  significantly  different  from  those  measured  experimentally. 
This  was  probably  caused  by  a  reduction  of  the  thermal 
conductivity  of  the  composite  material  near  the  heat-transfer 
surface  due  to  both  the  wall  thermal  resistance  and  the  contact 
thermal  resistance  [190]. 

The  synthesis  for  porous  and  fibrous  materials  is  given  in 
Table  7. 


124 


Y.  Dutil  et  al.  / Renewable  and  Sustainable  Energy  Reviews  15  (201 1)  112-130 


Table  7 

Models  for  porous  and  fibrous  materials. 


Model 

Material 

Numerical 

formulation 

Comment 

Validation 

Erk  and  Dudukovic  [183] 

n-Octadecane 

MM 

ID 

Second  law  efficiency  calculated 

Yes 

Mesalhy  et  al.  [184,185] 

FD 

2D 

Carbon  foam  impregnated  with  PCM 

Khillarkar  et  al.  [187] 

FG 

Fukai  et  al.  [188,189] 

n-octadecane 

FD 

2D 

Carbon-fiber  brushes  in  PCM 

Yes 

FG 

Hamada  et  al.  [190] 

n-octadecane 

FD 

2D 

Carbon-fiber  chips  in  PCM 

Yes 

FG 

Elgafy  and  Lafdi  [191] 

Parafin  wax 

A 

ID 

Carbon  nanofibers  in  paraffin  wax 

Lafdi  and  Matzek  [192,191] 

Frusteri  et  al.  [193] 

37:26:38 

NH4NO3 

Mg(N03)2-6H20 

MgCl2-6H20 

2D 

Carbon  fibers  in  PCM 

Yes 

4.7.  Slurry 

Like  porous  and  fibrous  material,  slurries  increase  heat-transfer 
performance  and  capacity  of  fluid.  This  enhancement  improves  the 
overall  efficiency  of  cooling  or  heating  system.  Reader  should 
consult  the  review  of  Zhang  et  al.  [194]  to  learn  more  about 
properties  and  applications  of  phase-change  material  slurries. 

This  research  field  was  first  explored  by  Charunyakor  et  al. 
[195].  They  formulated  a  model  for  heat  transfer  of  microencap¬ 
sulated  phase-change  material  slurry  flow  in  circular  ducts.  Heat 
generation  or  absorption  due  to  phase  change  in  the  particles  was 
included  in  the  energy  equation  as  a  source  term.  The  enhance¬ 
ment  of  thermal  conductivity  due  to  the  particle/fluid  interactions 
was  also  taken  into  consideration.  Results  of  this  model  have 
demonstrated  that  Stefan  number  and  PCM  microcapsule  concen¬ 
tration  are  the  most  important  parameters  for  the  system 
performances. 

Zhang  and  Faghri  [196]  studied  laminar  forced  convection  heat 
transfer  of  a  microencapsulated  phase-change  material  in  a 
circular  tube  with  constant  heat  flux.  Melting  in  the  microcapsule 
was  solved  by  a  temperature  transforming  model  instead  of  a 
quasi-steady  model.  This  model  was  validated  with  results  of 
Charunyakorn  et  al.  [195]  and  Goel  et  al.  [197],  This  analysis 
pointed  out  that  the  most  significant  reason  for  the  large  difference 
between  Goel  et  al.  [197]  experimental  results  and  Charunyakorn 
et  al.  [195]  numerical  results  is  that  melting  takes  place  over  a 
range  of  temperatures  below  the  melting  point.  It  also  put  into 
evidence  the  importance  of  shell  heat-transfer  properties  to 
explain  this  discrepancy. 

Alisetti  et  al.  [198,199]  introduced  an  effective  heat  capacity 
model  for  heat  transfer  in  PCM  slurries.  Their  objective  was  to 
address  the  shortcoming  of  previous  models  using  complicated 
source  terms  or  special  analytical  techniques  or  variation  of  the 
specific  heat  of  material  with  temperature.  This  new  formulation 
was  easier  to  implement  in  standard  computer  fluid  dynamic 
packages. 

Roy  and  Avanic  [200]  developed  a  turbulent  heat-transfer 
model  based  on  the  effective  specific  heat  capacity  for  PCM  slurry 
in  a  circular  tube  with  constant  wall  heat  flux.  Their  model  was 
compared  with  previously  published  numerical  [237,238]  and 
experimental  data  [239-241  ].  It  demonstrated  that  the  bulk  Stefan 
number,  the  melt  temperature  range,  and  the  degree  of  subcooling 
are  the  driving  parameters  of  the  system. 

Hu  and  Zhang  [201  ]  presented  an  analysis  of  the  forced 
convective  heat-transfer  enhancement  of  microencapsulated 
phase-change  material  slurries  flowing  through  a  circular  tube 
with  constant  heat  flux.  The  model  used  an  effective  specific  heat 
capacity  and  was  validated  with  the  results  of  Goel  et  al.  [197].  It 
was  found  that  the  conventional  Nusselt  number  correlations  for 
internal  flow  of  single  phase  fluids  are  not  suitable  to  accurately 


describe  the  heat-transfer  enhancement  with  microencapsulated 
phase-change  material  suspensions,  and  a  modification  is  intro¬ 
duced  that  enables  evaluation  for  the  convective  heat  transfer  of 
internal  flows.  Results  from  the  model  show  that  Stephan  number 
and  volumetric  concentration  of  microcapsules  are  the  most 
important  parameters  governing  the  heat-transfer  enhancement 
of  phase-change  slurries.  However,  initial  supercooling,  phase- 
change  temperature  range  and  microcapsule  diameter  also 
influence  the  degree  of  heat-transfer  enhancement.  The  enhance¬ 
ment  increases  with  decreasing  initial  supercooling  and/or  phase- 
change  temperature  range,  and  with  microcapsule  diameter. 
Therefore,  the  degree  of  enhancement  in  a  small-diameter  tube 
may  be  better  than  in  a  large  pipe  for  a  given  phase-change  slurry. 

Hao  and  Tao  [202]  developed  a  model  for  simulation  of  the 
laminar  hydrodynamic  and  heat-transfer  characteristics  of  sus¬ 
pension  flow  with  micro-nano-size  PCM  particles  in  a  micro- 
channel.  The  model  solves  simultaneously  the  two-phase  equa¬ 
tions  of  mass,  momentum,  and  energy.  This  has  allowed  the 
investigation  of  heat-transfer  enhancement  limits,  due  to  non- 
thermal  equilibrium  melting  of  PCM,  existence  of  a  particle- 
depleted  layer,  particle-wall  interaction,  phase-change  region 
propagation,  and  local  heat-transfer  coefficients  for  both  constant 
wall  heat  flux  and  constant  wall  temperature  boundary  conditions. 
The  same  group  [203]  developed  a  two-phase,  non  thermal 
equilibrium-based  model  for  the  numerical  simulation  of  laminar 
flow  and  heat-transfer  characteristics  of  PCM  slurry  in  a  micro- 
channel.  The  model  was  validated  against  the  experimental  results 
of  Goel  et  al.  [197].  The  results  show  that  for  a  given  Reynolds 
number,  there  exists  an  optimal  heat  flux  under  which  the 
effectiveness  is  the  highest. 

Ho  et  al.  [204]  examined  the  effects  of  wall  conduction  on  the 
heat-transfer  characteristics  of  solid-liquid  phase-change  material 
suspension  flow.  A  mixture  continuum  approach  was  adopted  in 
the  formulation  of  the  energy  equation,  with  an  enthalpy  model 
describing  the  phase-change  process  in  the  phase-change  material 
particles.  From  numerical  simulations  via  the  finite-volume 
approach,  it  was  found  that  the  conduction  heat-transfer 
propagating  along  the  pipe  wall  resulted  in  significant  preheating 
of  the  suspension  flow  in  the  region  upstream  of  the  heated 
section,  where  melting  of  the  particles  may  occur  and  therefore  the 
contribution  of  the  latent  heat  transfer  to  convection  heat 
dissipation  over  the  heated  section  was  markedly  attenuated. 

Inaba  et  al.  [205]  studied  the  fluid  flow  and  heat-transfer 
characteristics  for  Rayleigh-Benard  natural  convection  of  non- 
Newtonian  phase-change  material.  From  their  model,  they  found 
that  Rayleigh  number,  Prandtl  number  and  aspect  ratio  were  the 
most  important  factors  for  evaluating  the  natural  convection  in 
enclosures  for  most  of  Newtonian  and  non-Newtonian  fluids. 
However,  some  modifications  are  necessary  for  evaluating  the 
natural  convection  in  a  PCM  slurry.  In  consequence,  a  modified 
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Model 

Material 

Numerical 

formulation 

Comment 

Validation 

Charunyakor  et  al.  [195] 

FD 

ID 

Circular  ducts 

Ahuja  [236] 

Zhang  and  Faghri  [196] 

FD  FG 

ID 

Circular  tubes,  constant  heat  flux 

Experimental 

Goel  et  al.  [197] 

Alisetti  et  al.  [198,199] 

CFD 

Effective  heat  capacity 

Roy  and  Avanic  [200] 

10%  n-hexadecane 

FD 

3D 

Turbulent  heat  transfer 

Numerical  Petukhov  [237] 

in  water 

FG 

Sparrow  [238] 

Experimental 

Hartnett  [239] 

Choi  [240,241] 

Hu  et  al.  [201,226] 

FD 

2D 

Circular  tubes,  constant  heat  flux 

Experimental 

FG 

Goel  et  al.  [197] 

Ho  et  al.  [204] 

FD 

Micro  and  nano  particles 

Experimental  [197] 

FG 

Numerical  [195,196] 

Inaba  et  al.  [205,206] 

30:5:65 

FD 

2D 

Natural  convection,  Newtonian  and 

Yes 

Parafin  surfactant 

FG 

non-Newtonian  fluids 

Numerical 

water 

Vahl  Davis  and  Jones  [242] 
Ozoe  and  Churchill  [243] 

Cassidy  [207] 

Octacosane 

FD 

2D 

Mixed  convection 

Yes 

FG 

Royon  and  Guiffant  [208] 

FD 

ID 

Millimetric  particles  in  oil 

Numerical  Bird  et  al.  [244] 

FG 

Sabbah  et  al.  [209] 

CFD 

3D 

3D  single  phase  flow,  FLUENT 

Experimental 

Goel  et  al.  [197] 

Kuravi  et  al.  [210] 

FE 

3D 

Bulk  fluid  with  adaptative  specific  heat 

Experimental 

Goel  et  al.  [197] 

Zhang  et  al.  [211,212] 

FD 

2D 

Analogy  with  thermal  sources 

Charunyakorn  et  al.  [195] 

FG 

Alisetti  and  Roy  [199] 

Sun  et  al.  [213] 

Paraffin  wax 

2D 

Neutrally  buoyant  ceramics 

Experimental 

Jones  et  al.  [118] 

Stefan  number  was  defined  to  address  this  issue.  The  same  team 
[206]  produced  a  two-dimensional  numerical  simulation  of 
natural  convection  in  a  rectangular  enclosure  with  non-Newtonian 
PCM  microcapsulate  slurry.  The  formulation  of  the  numerical 
model  has  been  done  using  the  finite-volume  method.  This  study 
has  demonstrated  that  heat-transfer  coefficients  of  the  PCM  slurry 
with  phase  change  are  higher  than  those  without  phase  change.  It 
was  found  that  the  natural  convection  effect  is  intensified,  and  the 
heat  transfer  is  enhanced  in  the  case  of  the  PCM  slurry,  due  to  the 
contribution  of  the  latent  heat  transfer. 

Cassidy  [207]  studied  numerically  and  experimentally  a  micro 
PCM  fluid  in  a  circular  tube  with  mixed  convection.  The  PCM  was 
octacosane  encapsulated  by  a  polyethylene  shell  to  form  a 
spherical  particle  with  an  average  diameter  of  20  p,m.  The  micro 
PCM  particles  were  suspended  in  a  50/50  ethylene  glycol  water 
mixture.  Experimental  measurements  were  made  of  the  tube  outer 
wall  at  the  top  and  bottom  of  the  copper  tube.  Numerically  an 
incompressible  flow  model  was  used  with  an  Eulerian-Eulerian 
method  to  solve  the  two-phase  momentum  and  energy  equations. 
The  numerical  model  was  validated  using  experimental  data  of 
single  phase  flow  with  mixed  convection  available  in  the  literature 
and  was  also  verified  and  thermal  results  of  both  single  phase  and 
two-phase  flow  from  the  experimental  work  of  the  same  author. 

Royon  and  Cuiffant  [208]  constructed  a  model  describing  the 
thermal  behavior  of  slurry  of  phase-change  material  flow  with 
millimetric  particles  dispersed  in  oil  in  a  circular  duct.  They 
formulated  a  phenomenological  equation  which  takes  into  account 
the  heat  generation  due  to  phase  change  in  the  particles.  A  plateau 
of  the  temperature  appeared  in  all  of  the  simulations  where  the 
PCM  particles  were  considered.  The  length  of  this  plateau  is 
strongly  dependent  on  both  the  weight  fraction  of  the  particles  and 
the  duct  wall  temperature. 

Sabbah  et  al.  [209]  has  investigated  the  influence  of  micro- 
encapsulated  phase-change  material  on  the  thermal  and  hydraulic 
performance  of  micro-channel  heat  sinks.  They  developed  a  three¬ 


dimensional,  one-phase,  laminar  flow  model  of  a  rectangular 
channel  using  water  slurry  of  MEPCM  with  temperature  depen¬ 
dent  physical  properties.  The  conservation  equation  was  solved 
using  FLUENT  6.2.  The  numerical  results  have  been  compared  to 
experimental  data  of  Goel  et  al.  [197]  once  adapted  for  the  circular 
geometry  of  this  experiment.  The  model  results  showed  a 
significant  increase  in  the  heat-transfer  coefficient  that  is  mainly 
dependant  on  the  channel  inlet  and  outlet  temperatures  and  the 
selected  PCM  melting  temperature.  The  enhancement  is  higher  for 
low  concentration  of  PCM  in  slurry  due  to  the  increase  of  slurry 
viscosity  with  its  concentration. 

Kuravi  et  al.  [210]  produced  a  numerical  investigation  of  the 
performance  of  a  nano-encapsulated  phase-change  material  slurry 
in  a  manifold  micro-channel  heat  sink.  The  slurry  was  modeled  as  a 
bulk  fluid  with  varying  specific  heat.  The  temperature  field  inside 
the  channel  wall  is  solved  three  dimensionally  and  is  coupled  with 
the  three-dimensional  velocity  and  temperature  fields  of  the  fluid. 

Zhang  et  al.  [211]  analyzed  the  convective  heat-transfer 
enhancement  mechanism  of  microencapsulated  PCM  slurries 
based  on  the  analogy  between  convective  heat  transfer  and 
thermal  conduction  with  thermal  sources.  The  heat-transfer 
enhancement  for  laminar  flow  in  a  circular  tube  with  constant 
wall  temperature  is  analyzed  using  an  effective  specific  heat 
capacity  model.  This  model  was  validated  by  comparison  with  the 
result  of  Charunyakorn  et  al.  [195]  and  Alisetti  and  Roy  [199].  The 
numerical  simulation  results  have  pointed  out  the  Stephan 
number  and  the  phase-change  temperature  range  as  the  most 
important  parameters  influencing  the  Nusselt  number  fluctuation 
profile  of  phase-change  slurry.  The  fluctuation  amplitude  increases 
with  decreasing  phase-change  temperature  range.  The  fluctuation 
amplitude  and  range  increase  with  decreasing  Stephan  number 
when  phase-change  happens.  Later,  members  of  the  same  team 
[212]  studied  experimentally  and  numerically  characteristics  of 
microencapsulated  PCM  flowing  in  a  circular  tube.  The  enhanced 
convective  heat-transfer  mechanism  was  analyzed  by  using  the 
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enthalpy  model.  Three  kinds  of  fluid-pure  water,  micro-particle 
slurry  and  MPCMs  were  numerically  investigated.  They  had 
essentially  reached  the  same  conclusions. 

Using  the  same  setup  as  Jones  et  al.  [119],  Sun  et  al.  [213] 
studied  the  behavior  of  slurry  composed  of  neutrally  buoyant 
ceramic  hollow  spheres  suspended  in  a  paraffin  wax.  The 
numerical  analysis  employed  a  particle-diffusive  model  and  the 
enthalpy  method.  Reasonable  agreement  was  obtained  between 
the  experiments  and  numerical  predictions.  Nevertheless,  flow  and 
heat-transfer  characteristics  were  greatly  altered  due  to  the 
presence  of  the  solid  particles  and  that  the  particle-diffusive 
model  was  insufficient  to  describe  the  particle  migration  during 
melting. 

A  synthesis  of  the  numerical  study  of  slurries  is  provided  in 
Table  8. 

5.  Conclusion 

5.3.  Content 

This  paper  presented  an  exhaustive  review  of  numerical 
methods  applied  to  the  solutions  of  heat-transfer  problems 
involving  phase-change  materials  for  thermal  energy  storage. 
The  review  is  a  model  collection  of  fundamental  and  most  recent 
works  published  on  the  subject.  This  survey  is  organized  according 
to  the  problem  geometry  (Cartesian,  spherical,  and  cylindrical)  and 
specific  configurations  or  applications  (packed  beds,  finned 
surfaces,  porous  and  fibrous  materials,  slurries).  The  authors  do 
not  claim  anything  about  the  completeness  of  the  review  as  some 
papers  may  have  been  unfortunately  neglected.  The  authors 
apology  for  any  crucial  omission  and  welcome  all  authors  to 
contact  them  to  complete  the  review. 

5.2.  Analysis:  numerical  issues  and  validation 

As  a  general  conclusion,  one  can  readily  observe  that  the  older 
studies  all  had  experimental  counterparts  to  validate  the  modeling 
of  the  problem  with  an  appropriate  set  of  equations  and  the 
subsequent  formulation  of  a  numerical  method  to  solve  the 
relevant  sets  of  discretized  equations.  Many  of  these  early  studies 
also  involved  analytical  solutions  used  to  validate  the  model  for 
selected  problems  that  admit  closed-form  solutions. 

As  time  passes,  the  authors  rely  more  and  more  on  other 
studies,  mostly  other  numerical  studies,  to  validate  their  own 
numerical  results.  And  it  is  somewhat  interesting  to  mention  that 
among  the  250+  references  cited  here,  in  only  one  the  authors 
stated  that  the  results  were  not  "in  good  agreement  with  those 
found  in  the  references”.  Many  of  the  recent  studies  discuss  their 
results  qualitatively  only  as  the  comparison  with  a  graph  taken 
from  a  publication  may  be  somehow  hazardous. 

In  recent  studies,  the  proportion  of  analyses  which  rely  on 
commercial  codes  increases  and  the  discussions  that  pertain  to 
stability,  convergence,  grid  independence  and  other  related 
numerical  issues  diminishes. 

On  the  other  hand,  we  found  that  some  laboratories  are  involved 
in  LHTES  for  decades  and  in  these  labs  there  is  a  mature  culture  of 
experimental,  analytical,  and  numerical  methods  for  the  analysis. 

5.3.  Analysis:  type  of  materials  and  configurations 

We  would  like  to  warn  readers  about  the  risk  involved  when 
extrapolating  from  numerical  and  experimental  studies  and 
making  general  conclusion  from  them.  For  example,  convection 
is  strongly  sensitive  to  the  geometry  and  the  size  of  the  enclosure 
as  to  the  viscosity  and  thermal  conductivity  of  PCM.  Thermal 
conductivity  of  the  enclosure  and  the  PCM  solid-liquid  density  and 


conductivity  contrasts  also  influence  the  thermal  behavior  in 
important  and  counter-intuitive  ways. 

There  is  an  incentive  to  mix  several  materials  with  different 
phase-change  temperature  to  allow  for  more  heat  recovery.  Packed 
beds,  finned  geometries,  porous  and  fibrous  material-based 
applications,  and  slurries  are  the  particular  configurations  for 
most  studies  that  pertain  to  systems  involving  PCMs.  The  latter 
receives  more  attention  these  days  as  a  result  of  the  incentive  to 
study  nano-heat  transfer  in  OECD  countries. 

However,  it  would  have  been  interesting  in  these  studies  to 
discuss  the  issue  of  pumping  power  or  the  effect  of  the  micro-  or 
nano-capsules  of  PCMs  on  the  apparent  or  effective  viscosity  of  the 
fluids.  None  of  the  authors  mention  this  critical  issue  when  it 
comes  to  establish  the  net  energy  balance  of  operation. 

Moreover,  in  several  studies,  carbon  fibers  are  used  to  enhance 
conduction  in  PCMs  (it  is  a  well  known  fact  that  their  average 
conductivity  is  poor)  but  these  studies  are  silent  about  the  energy 
cost  of  production  of  those  fibers. 

Another  important  aspect  often  neglected  when  modeling 
LHTES  is  the  impact  of  thermal  dilatation  of  PCM,  which  must  be 
taken  into  account  into  the  conception  of  experimental  and 
operational  systems.  In  both  case,  geometry  differs  from  the 
idealized  version  used  in  most  models,  with  poorly  documented 
impact  on  performances. 

5.4.  Analysis:  economic  impact 

This  issue  is  not  treated  at  all  in  any  of  the  paper.  The  use  of  a 
PCM-based  system  compared  to  a  standard  sensible  heat  system  is 
not  discussed.  The  whole  issue  of  storage  in  a  technico-economical 
context  is  not  found  in  the  reviewed  papers. 

5.5.  Concluding  remarks 

Models  are  now  established  for  most  applications  of  PCM 
materials.  In  many  cases,  simplified  models  or  analytic  expressions 
or  correlation  functions  have  been  developed  for  practical 
guidance  in  optimizing  design  of  systems  using  PCMs.  They 
should  be  used  as  guidelines  when  adjusting  the  implementation 
of  a  numerical  method  to  be  used  in  the  design  of  systems. 

The  validation  of  numerical  predictions  should  always  be 
performed  by  means  of  appropriate  experimental  data.  The 
experimental  data  should  in  turn  be  accompanied  by  a  suitable 
uncertainty  analysis  while  the  predictions  should  first  address  the 
issues  of  stability,  convergence,  etc.  One  should  also  keep  in  mind 
that  the  so-called  "governing”  equations  do  not  govern  but 
describe  the  physics  of  the  problem.  This  semantic  distinction 
should  avoid  comments  found  in  a  paper  into  which  the  authors 
indicate  that  experiments  are  not  in  agreement  with  the 
predictions  instead  of  the  opposite. 

Life  cycle  analysis,  both  economical  and  energetic,  should  be 
performed  on  systems  to  determine  into  which  conditions  energy 
storage  systems  involving  PCMs  should  be  used. 

Despite  these  criticisms,  the  actual  situation  which  finds 
climate  changes,  security  of  the  supplies,  and  fossil  fuels  depletion 
at  the  heart  of  the  economic  discussions  will  ensure  a  place  for 
PCMs  in  the  global  energy  efficiency  policies  of  the  world  to  come. 
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